

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">


<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    
    <title>agpy.mpfit_custom &mdash; agpy 0.1.2 documentation</title>
    
    <link rel="stylesheet" href="../../_static/extra.css" type="text/css" />
    <link rel="stylesheet" href="../../_static/pygments.css" type="text/css" />
    
    <script type="text/javascript">
      var DOCUMENTATION_OPTIONS = {
        URL_ROOT:    '../../',
        VERSION:     '0.1.2',
        COLLAPSE_INDEX: false,
        FILE_SUFFIX: '.html',
        HAS_SOURCE:  true
      };
    </script>
    <script type="text/javascript" src="../../_static/jquery.js"></script>
    <script type="text/javascript" src="../../_static/underscore.js"></script>
    <script type="text/javascript" src="../../_static/doctools.js"></script>
    <link rel="top" title="agpy 0.1.2 documentation" href="../../index.html" />
    <link rel="up" title="agpy" href="../agpy.html" />
     
    <script type="text/javascript">

      var _gaq = _gaq || [];
      _gaq.push(['_setDomainName', 'pyspeckit.bitbucket.org']);
      _gaq.push(['_setAllowHash', false]);
      _gaq.push(['_trackPageview']);


    </script>
    <link rel="stylesheet" type="text/css" href="../../_static/extra.css" />
  </head>
  <body>
    <div class="header-wrapper">
      <div class="header">
        <h1><a href="../../index.html">agpy 0.1.2 documentation</a></h1>
        <div class="rel">
          <a href="http://agpy.googlecode.com">agpy Home </a>  |
          <a href=../../index.html>Docs Home </a>  |
          <a href="http://code.google.com/p/agpy/w/list">Wiki</a>  |
          <a href=../../search.html>Search </a>
        </div>
       </div>
    </div>

    <div class="content-wrapper">
      <div class="content">
        <div class="document">
            
      <div class="documentwrapper">
        <div class="bodywrapper">
          <div class="body">
            
  <h1>Source code for agpy.mpfit_custom</h1><div class="highlight"><pre>
<span class="sd">&quot;&quot;&quot;</span>
<span class="sd">Perform Levenberg-Marquardt least-squares minimization, based on MINPACK-1.</span>

<span class="sd">                                   AUTHORS</span>
<span class="sd">  The original version of this software, called LMFIT, was written in FORTRAN</span>
<span class="sd">  as part of the MINPACK-1 package by XXX.</span>

<span class="sd">  Craig Markwardt converted the FORTRAN code to IDL.  The information for the</span>
<span class="sd">  IDL version is:</span>
<span class="sd">     Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770</span>
<span class="sd">     craigm@lheamail.gsfc.nasa.gov</span>
<span class="sd">     UPDATED VERSIONs can be found on my WEB PAGE:</span>
<span class="sd">        http://cow.physics.wisc.edu/~craigm/idl/idl.html</span>
<span class="sd">    </span>
<span class="sd">  Mark Rivers created this Python version from Craig&#39;s IDL version.</span>
<span class="sd">    Mark Rivers, University of Chicago</span>
<span class="sd">    Building 434A, Argonne National Laboratory</span>
<span class="sd">    9700 South Cass Avenue, Argonne, IL 60439</span>
<span class="sd">    rivers@cars.uchicago.edu</span>
<span class="sd">    Updated versions can be found at http://cars.uchicago.edu/software</span>
<span class="sd"> </span>
<span class="sd"> Sergey Koposov converted the Mark&#39;s Python version from Numeric to numpy</span>
<span class="sd">    Sergey Koposov, Max Planck Institute for Astronomy</span>
<span class="sd">    Heidelberg, Germany, D-69117</span>
<span class="sd">    koposov@mpia.de</span>
<span class="sd">    Updated versions can be found at http://code.google.com/p/astrolibpy/source/browse/trunk/</span>

<span class="sd">                                 DESCRIPTION</span>

<span class="sd"> MPFIT uses the Levenberg-Marquardt technique to solve the</span>
<span class="sd"> least-squares problem.  In its typical use, MPFIT will be used to</span>
<span class="sd"> fit a user-supplied function (the &quot;model&quot;) to user-supplied data</span>
<span class="sd"> points (the &quot;data&quot;) by adjusting a set of parameters.  MPFIT is</span>
<span class="sd"> based upon MINPACK-1 (LMDIF.F) by More&#39; and collaborators.</span>

<span class="sd"> For example, a researcher may think that a set of observed data</span>
<span class="sd"> points is best modelled with a Gaussian curve.  A Gaussian curve is</span>
<span class="sd"> parameterized by its mean, standard deviation and normalization.</span>
<span class="sd"> MPFIT will, within certain constraints, find the set of parameters</span>
<span class="sd"> which best fits the data.  The fit is &quot;best&quot; in the least-squares</span>
<span class="sd"> sense; that is, the sum of the weighted squared differences between</span>
<span class="sd"> the model and data is minimized.</span>

<span class="sd"> The Levenberg-Marquardt technique is a particular strategy for</span>
<span class="sd"> iteratively searching for the best fit.  This particular</span>
<span class="sd"> implementation is drawn from MINPACK-1 (see NETLIB), and is much faster</span>
<span class="sd"> and more accurate than the version provided in the Scientific Python package</span>
<span class="sd"> in Scientific.Functions.LeastSquares.</span>
<span class="sd"> This version allows upper and lower bounding constraints to be placed on each</span>
<span class="sd"> parameter, or the parameter can be held fixed.</span>

<span class="sd"> The user-supplied Python function should return an array of weighted</span>
<span class="sd"> deviations between model and data.  In a typical scientific problem</span>
<span class="sd"> the residuals should be weighted so that each deviate has a</span>
<span class="sd"> gaussian sigma of 1.0.  If X represents values of the independent</span>
<span class="sd"> variable, Y represents a measurement for each value of X, and ERR</span>
<span class="sd"> represents the error in the measurements, then the deviates could</span>
<span class="sd"> be calculated as follows:</span>

<span class="sd">   DEVIATES = (Y - F(X)) / ERR</span>

<span class="sd"> where F is the analytical function representing the model.  You are</span>
<span class="sd"> recommended to use the convenience functions MPFITFUN and</span>
<span class="sd"> MPFITEXPR, which are driver functions that calculate the deviates</span>
<span class="sd"> for you.  If ERR are the 1-sigma uncertainties in Y, then</span>

<span class="sd">   TOTAL( DEVIATES^2 )</span>

<span class="sd"> will be the total chi-squared value.  MPFIT will minimize the</span>
<span class="sd"> chi-square value.  The values of X, Y and ERR are passed through</span>
<span class="sd"> MPFIT to the user-supplied function via the FUNCTKW keyword.</span>

<span class="sd"> Simple constraints can be placed on parameter values by using the</span>
<span class="sd"> PARINFO keyword to MPFIT.  See below for a description of this</span>
<span class="sd"> keyword.</span>

<span class="sd"> MPFIT does not perform more general optimization tasks.  See TNMIN</span>
<span class="sd"> instead.  MPFIT is customized, based on MINPACK-1, to the</span>
<span class="sd"> least-squares minimization problem.</span>


<span class="sd">                               USER FUNCTION</span>

<span class="sd"> The user must define a function which returns the appropriate</span>
<span class="sd"> values as specified above.  The function should return the weighted</span>
<span class="sd"> deviations between the model and the data.  It should also return a status</span>
<span class="sd"> flag and an optional partial derivative array.  For applications which</span>
<span class="sd"> use finite-difference derivatives -- the default -- the user</span>
<span class="sd"> function should be declared in the following way:</span>

<span class="sd">   def myfunct(p, fjac=None, x=None, y=None, err=None)</span>
<span class="sd">    # Parameter values are passed in &quot;p&quot;</span>
<span class="sd">    # If fjac==None then partial derivatives should not be</span>
<span class="sd">    # computed.  It will always be None if MPFIT is called with default</span>
<span class="sd">    # flag.</span>
<span class="sd">    model = F(x, p)</span>
<span class="sd">    # Non-negative status value means MPFIT should continue, negative means</span>
<span class="sd">    # stop the calculation.</span>
<span class="sd">    status = 0</span>
<span class="sd">    return([status, (y-model)/err]</span>

<span class="sd"> See below for applications with analytical derivatives.</span>

<span class="sd"> The keyword parameters X, Y, and ERR in the example above are</span>
<span class="sd"> suggestive but not required.  Any parameters can be passed to</span>
<span class="sd"> MYFUNCT by using the functkw keyword to MPFIT.  Use MPFITFUN and</span>
<span class="sd"> MPFITEXPR if you need ideas on how to do that.  The function *must*</span>
<span class="sd"> accept a parameter list, P.</span>

<span class="sd"> In general there are no restrictions on the number of dimensions in</span>
<span class="sd"> X, Y or ERR.  However the deviates *must* be returned in a</span>
<span class="sd"> one-dimensional Numeric array of type Float.</span>

<span class="sd"> User functions may also indicate a fatal error condition using the</span>
<span class="sd"> status return described above. If status is set to a number between</span>
<span class="sd"> -15 and -1 then MPFIT will stop the calculation and return to the caller.</span>


<span class="sd">                            ANALYTIC DERIVATIVES</span>

<span class="sd"> In the search for the best-fit solution, MPFIT by default</span>
<span class="sd"> calculates derivatives numerically via a finite difference</span>
<span class="sd"> approximation.  The user-supplied function need not calculate the</span>
<span class="sd"> derivatives explicitly.  However, if you desire to compute them</span>
<span class="sd"> analytically, then the AUTODERIVATIVE=0 keyword must be passed to MPFIT.</span>
<span class="sd"> As a practical matter, it is often sufficient and even faster to allow</span>
<span class="sd"> MPFIT to calculate the derivatives numerically, and so</span>
<span class="sd"> AUTODERIVATIVE=0 is not necessary.</span>

<span class="sd"> If AUTODERIVATIVE=0 is used then the user function must check the parameter</span>
<span class="sd"> FJAC, and if FJAC!=None then return the partial derivative array in the</span>
<span class="sd"> return list.</span>
<span class="sd">   def myfunct(p, fjac=None, x=None, y=None, err=None)</span>
<span class="sd">    # Parameter values are passed in &quot;p&quot;</span>
<span class="sd">    # If FJAC!=None then partial derivatives must be comptuer.</span>
<span class="sd">    # FJAC contains an array of len(p), where each entry</span>
<span class="sd">    # is 1 if that parameter is free and 0 if it is fixed.</span>
<span class="sd">    model = F(x, p)</span>
<span class="sd">    Non-negative status value means MPFIT should continue, negative means</span>
<span class="sd">    # stop the calculation.</span>
<span class="sd">    status = 0</span>
<span class="sd">    if (dojac):</span>
<span class="sd">       pderiv = zeros([len(x), len(p)], Float)</span>
<span class="sd">       for j in range(len(p)):</span>
<span class="sd">         pderiv[:,j] = FGRAD(x, p, j)</span>
<span class="sd">    else:</span>
<span class="sd">       pderiv = None</span>
<span class="sd">    return([status, (y-model)/err, pderiv]</span>

<span class="sd"> where FGRAD(x, p, i) is a user function which must compute the</span>
<span class="sd"> derivative of the model with respect to parameter P[i] at X.  When</span>
<span class="sd"> finite differencing is used for computing derivatives (ie, when</span>
<span class="sd"> AUTODERIVATIVE=1), or when MPFIT needs only the errors but not the</span>
<span class="sd"> derivatives the parameter FJAC=None.</span>

<span class="sd"> Derivatives should be returned in the PDERIV array. PDERIV should be an m x</span>
<span class="sd"> n array, where m is the number of data points and n is the number</span>
<span class="sd"> of parameters.  dp[i,j] is the derivative at the ith point with</span>
<span class="sd"> respect to the jth parameter.</span>

<span class="sd"> The derivatives with respect to fixed parameters are ignored; zero</span>
<span class="sd"> is an appropriate value to insert for those derivatives.  Upon</span>
<span class="sd"> input to the user function, FJAC is set to a vector with the same</span>
<span class="sd"> length as P, with a value of 1 for a parameter which is free, and a</span>
<span class="sd"> value of zero for a parameter which is fixed (and hence no</span>
<span class="sd"> derivative needs to be calculated).</span>

<span class="sd"> If the data is higher than one dimensional, then the *last*</span>
<span class="sd"> dimension should be the parameter dimension.  Example: fitting a</span>
<span class="sd"> 50x50 image, &quot;dp&quot; should be 50x50xNPAR.</span>


<span class="sd">           CONSTRAINING PARAMETER VALUES WITH THE PARINFO KEYWORD</span>

<span class="sd"> The behavior of MPFIT can be modified with respect to each</span>
<span class="sd"> parameter to be fitted.  A parameter value can be fixed; simple</span>
<span class="sd"> boundary constraints can be imposed; limitations on the parameter</span>
<span class="sd"> changes can be imposed; properties of the automatic derivative can</span>
<span class="sd"> be modified; and parameters can be tied to one another.</span>

<span class="sd"> These properties are governed by the PARINFO structure, which is</span>
<span class="sd"> passed as a keyword parameter to MPFIT.</span>

<span class="sd"> PARINFO should be a list of dictionaries, one list entry for each parameter.</span>
<span class="sd"> Each parameter is associated with one element of the array, in</span>
<span class="sd"> numerical order.  The dictionary can have the following keys</span>
<span class="sd"> (none are required, keys are case insensitive):</span>

<span class="sd">    &#39;value&#39; - the starting parameter value (but see the START_PARAMS</span>
<span class="sd">             parameter for more information).</span>

<span class="sd">    &#39;fixed&#39; - a boolean value, whether the parameter is to be held</span>
<span class="sd">             fixed or not.  Fixed parameters are not varied by</span>
<span class="sd">             MPFIT, but are passed on to MYFUNCT for evaluation.</span>

<span class="sd">    &#39;limited&#39; - a two-element boolean array.  If the first/second</span>
<span class="sd">               element is set, then the parameter is bounded on the</span>
<span class="sd">               lower/upper side.  A parameter can be bounded on both</span>
<span class="sd">               sides.  Both LIMITED and LIMITS must be given</span>
<span class="sd">               together.</span>

<span class="sd">    &#39;limits&#39; - a two-element float array.  Gives the</span>
<span class="sd">              parameter limits on the lower and upper sides,</span>
<span class="sd">              respectively.  Zero, one or two of these values can be</span>
<span class="sd">              set, depending on the values of LIMITED.  Both LIMITED</span>
<span class="sd">              and LIMITS must be given together.</span>

<span class="sd">    &#39;parname&#39; - a string, giving the name of the parameter.  The</span>
<span class="sd">               fitting code of MPFIT does not use this tag in any</span>
<span class="sd">               way.  However, the default iterfunct will print the</span>
<span class="sd">               parameter name if available.</span>

<span class="sd">    &#39;step&#39; - the step size to be used in calculating the numerical</span>
<span class="sd">            derivatives.  If set to zero, then the step size is</span>
<span class="sd">            computed automatically.  Ignored when AUTODERIVATIVE=0.</span>

<span class="sd">    &#39;mpside&#39; - the sidedness of the finite difference when computing</span>
<span class="sd">              numerical derivatives.  This field can take four</span>
<span class="sd">              values:</span>

<span class="sd">                 0 - one-sided derivative computed automatically</span>
<span class="sd">                 1 - one-sided derivative (f(x+h) - f(x)  )/h</span>
<span class="sd">                -1 - one-sided derivative (f(x)   - f(x-h))/h</span>
<span class="sd">                 2 - two-sided derivative (f(x+h) - f(x-h))/(2*h)</span>

<span class="sd">             Where H is the STEP parameter described above.  The</span>
<span class="sd">             &quot;automatic&quot; one-sided derivative method will chose a</span>
<span class="sd">             direction for the finite difference which does not</span>
<span class="sd">             violate any constraints.  The other methods do not</span>
<span class="sd">             perform this check.  The two-sided method is in</span>
<span class="sd">             principle more precise, but requires twice as many</span>
<span class="sd">             function evaluations.  Default: 0.</span>

<span class="sd">    &#39;mpmaxstep&#39; - the maximum change to be made in the parameter</span>
<span class="sd">                 value.  During the fitting process, the parameter</span>
<span class="sd">                 will never be changed by more than this value in</span>
<span class="sd">                 one iteration.</span>

<span class="sd">                 A value of 0 indicates no maximum.  Default: 0.</span>

<span class="sd">    &#39;tied&#39; - a string expression which &quot;ties&quot; the parameter to other</span>
<span class="sd">            free or fixed parameters.  Any expression involving</span>
<span class="sd">            constants and the parameter array P are permitted.</span>
<span class="sd">            Example: if parameter 2 is always to be twice parameter</span>
<span class="sd">            1 then use the following: parinfo(2).tied = &#39;2 * p(1)&#39;.</span>
<span class="sd">            Since they are totally constrained, tied parameters are</span>
<span class="sd">            considered to be fixed; no errors are computed for them.</span>
<span class="sd">            [ NOTE: the PARNAME can&#39;t be used in expressions. ]</span>

<span class="sd">    &#39;mpprint&#39; - if set to 1, then the default iterfunct will print the</span>
<span class="sd">               parameter value.  If set to 0, the parameter value</span>
<span class="sd">               will not be printed.  This tag can be used to</span>
<span class="sd">               selectively print only a few parameter values out of</span>
<span class="sd">               many.  Default: 1 (all parameters printed)</span>


<span class="sd"> Future modifications to the PARINFO structure, if any, will involve</span>
<span class="sd"> adding dictionary tags beginning with the two letters &quot;MP&quot;.</span>
<span class="sd"> Therefore programmers are urged to avoid using tags starting with</span>
<span class="sd"> the same letters; otherwise they are free to include their own</span>
<span class="sd"> fields within the PARINFO structure, and they will be ignored.</span>

<span class="sd"> PARINFO Example:</span>
<span class="sd"> parinfo = [{&#39;value&#39;:0., &#39;fixed&#39;:0, &#39;limited&#39;:[0,0], &#39;limits&#39;:[0.,0.]} </span>
<span class="sd">                                                for i in range(5)]</span>
<span class="sd"> parinfo[0][&#39;fixed&#39;] = 1</span>
<span class="sd"> parinfo[4][&#39;limited&#39;][0] = 1</span>
<span class="sd"> parinfo[4][&#39;limits&#39;][0]  = 50.</span>
<span class="sd"> values = [5.7, 2.2, 500., 1.5, 2000.]</span>
<span class="sd"> for i in range(5): parinfo[i][&#39;value&#39;]=values[i]</span>

<span class="sd"> A total of 5 parameters, with starting values of 5.7,</span>
<span class="sd"> 2.2, 500, 1.5, and 2000 are given.  The first parameter</span>
<span class="sd"> is fixed at a value of 5.7, and the last parameter is</span>
<span class="sd"> constrained to be above 50.</span>


<span class="sd">                                   EXAMPLE</span>

<span class="sd">   import mpfit</span>
<span class="sd">   import numpy.oldnumeric as Numeric</span>
<span class="sd">   x = arange(100, float)</span>
<span class="sd">   p0 = [5.7, 2.2, 500., 1.5, 2000.]</span>
<span class="sd">   y = ( p[0] + p[1]*[x] + p[2]*[x**2] + p[3]*sqrt(x) +</span>
<span class="sd">         p[4]*log(x))</span>
<span class="sd">   fa = {&#39;x&#39;:x, &#39;y&#39;:y, &#39;err&#39;:err}</span>
<span class="sd">   m = mpfit(&#39;myfunct&#39;, p0, functkw=fa)</span>
<span class="sd">   print &#39;status = &#39;, m.status</span>
<span class="sd">   if (m.status &lt;= 0): print &#39;error message = &#39;, m.errmsg</span>
<span class="sd">   print &#39;parameters = &#39;, m.params</span>

<span class="sd">   Minimizes sum of squares of MYFUNCT.  MYFUNCT is called with the X,</span>
<span class="sd">   Y, and ERR keyword parameters that are given by FUNCTKW.  The</span>
<span class="sd">   results can be obtained from the returned object m.</span>


<span class="sd">                            THEORY OF OPERATION</span>

<span class="sd">   There are many specific strategies for function minimization.  One</span>
<span class="sd">   very popular technique is to use function gradient information to</span>
<span class="sd">   realize the local structure of the function.  Near a local minimum</span>
<span class="sd">   the function value can be taylor expanded about x0 as follows:</span>

<span class="sd">      f(x) = f(x0) + f&#39;(x0) . (x-x0) + (1/2) (x-x0) . f&#39;&#39;(x0) . (x-x0)</span>
<span class="sd">             -----   ---------------   -------------------------------  (1)</span>
<span class="sd">     Order  0th       1st                     2nd</span>

<span class="sd">   Here f&#39;(x) is the gradient vector of f at x, and f&#39;&#39;(x) is the</span>
<span class="sd">   Hessian matrix of second derivatives of f at x.  The vector x is</span>
<span class="sd">   the set of function parameters, not the measured data vector.  One</span>
<span class="sd">   can find the minimum of f, f(xm) using Newton&#39;s method, and</span>
<span class="sd">   arrives at the following linear equation:</span>

<span class="sd">      f&#39;&#39;(x0) . (xm-x0) = - f&#39;(x0)                          (2)</span>

<span class="sd">   If an inverse can be found for f&#39;&#39;(x0) then one can solve for</span>
<span class="sd">   (xm-x0), the step vector from the current position x0 to the new</span>
<span class="sd">   projected minimum.  Here the problem has been linearized (ie, the</span>
<span class="sd">   gradient information is known to first order).  f&#39;&#39;(x0) is</span>
<span class="sd">   symmetric n x n matrix, and should be positive definite.</span>

<span class="sd">   The Levenberg - Marquardt technique is a variation on this theme.</span>
<span class="sd">   It adds an additional diagonal term to the equation which may aid the</span>
<span class="sd">   convergence properties:</span>

<span class="sd">      (f&#39;&#39;(x0) + nu I) . (xm-x0) = -f&#39;(x0)                (2a)</span>

<span class="sd">   where I is the identity matrix.  When nu is large, the overall</span>
<span class="sd">   matrix is diagonally dominant, and the iterations follow steepest</span>
<span class="sd">   descent.  When nu is small, the iterations are quadratically</span>
<span class="sd">   convergent.</span>

<span class="sd">   In principle, if f&#39;&#39;(x0) and f&#39;(x0) are known then xm-x0 can be</span>
<span class="sd">   determined.  However the Hessian matrix is often difficult or</span>
<span class="sd">   impossible to compute.  The gradient f&#39;(x0) may be easier to</span>
<span class="sd">   compute, if even by finite difference techniques.  So-called</span>
<span class="sd">   quasi-Newton techniques attempt to successively estimate f&#39;&#39;(x0)</span>
<span class="sd">   by building up gradient information as the iterations proceed.</span>

<span class="sd">   In the least squares problem there are further simplifications</span>
<span class="sd">   which assist in solving eqn (2).  The function to be minimized is</span>
<span class="sd">   a sum of squares:</span>

<span class="sd">       f = Sum(hi^2)                                         (3)</span>

<span class="sd">   where hi is the ith residual out of m residuals as described</span>
<span class="sd">   above.  This can be substituted back into eqn (2) after computing</span>
<span class="sd">   the derivatives:</span>

<span class="sd">       f&#39;  = 2 Sum(hi  hi&#39;)</span>
<span class="sd">       f&#39;&#39; = 2 Sum(hi&#39; hj&#39;) + 2 Sum(hi hi&#39;&#39;)                (4)</span>

<span class="sd">   If one assumes that the parameters are already close enough to a</span>
<span class="sd">   minimum, then one typically finds that the second term in f&#39;&#39; is</span>
<span class="sd">   negligible [or, in any case, is too difficult to compute].  Thus,</span>
<span class="sd">   equation (2) can be solved, at least approximately, using only</span>
<span class="sd">   gradient information.</span>

<span class="sd">   In matrix notation, the combination of eqns (2) and (4) becomes:</span>

<span class="sd">        hT&#39; . h&#39; . dx = - hT&#39; . h                         (5)</span>

<span class="sd">   Where h is the residual vector (length m), hT is its transpose, h&#39;</span>
<span class="sd">   is the Jacobian matrix (dimensions n x m), and dx is (xm-x0).  The</span>
<span class="sd">   user function supplies the residual vector h, and in some cases h&#39;</span>
<span class="sd">   when it is not found by finite differences (see MPFIT_FDJAC2,</span>
<span class="sd">   which finds h and hT&#39;).  Even if dx is not the best absolute step</span>
<span class="sd">   to take, it does provide a good estimate of the best *direction*,</span>
<span class="sd">   so often a line minimization will occur along the dx vector</span>
<span class="sd">   direction.</span>

<span class="sd">   The method of solution employed by MINPACK is to form the Q . R</span>
<span class="sd">   factorization of h&#39;, where Q is an orthogonal matrix such that QT .</span>
<span class="sd">   Q = I, and R is upper right triangular.  Using h&#39; = Q . R and the</span>
<span class="sd">   ortogonality of Q, eqn (5) becomes</span>

<span class="sd">        (RT . QT) . (Q . R) . dx = - (RT . QT) . h</span>
<span class="sd">                     RT . R . dx = - RT . QT . h         (6)</span>
<span class="sd">                          R . dx = - QT . h</span>

<span class="sd">   where the last statement follows because R is upper triangular.</span>
<span class="sd">   Here, R, QT and h are known so this is a matter of solving for dx.</span>
<span class="sd">   The routine MPFIT_QRFAC provides the QR factorization of h, with</span>
<span class="sd">   pivoting, and MPFIT_QRSOLV provides the solution for dx.</span>


<span class="sd">                                 REFERENCES</span>

<span class="sd">   MINPACK-1, Jorge More&#39;, available from netlib (www.netlib.org).</span>
<span class="sd">   &quot;Optimization Software Guide,&quot; Jorge More&#39; and Stephen Wright,</span>
<span class="sd">     SIAM, *Frontiers in Applied Mathematics*, Number 14.</span>
<span class="sd">   More&#39;, Jorge J., &quot;The Levenberg-Marquardt Algorithm:</span>
<span class="sd">     Implementation and Theory,&quot; in *Numerical Analysis*, ed. Watson,</span>
<span class="sd">     G. A., Lecture Notes in Mathematics 630, Springer-Verlag, 1977.</span>


<span class="sd">                           MODIFICATION HISTORY</span>

<span class="sd">   Translated from MINPACK-1 in FORTRAN, Apr-Jul 1998, CM</span>
<span class="sd"> Copyright (C) 1997-2002, Craig Markwardt</span>
<span class="sd"> This software is provided as is without any warranty whatsoever.</span>
<span class="sd"> Permission to use, copy, modify, and distribute modified or</span>
<span class="sd"> unmodified copies is granted, provided this copyright and disclaimer</span>
<span class="sd"> are included unchanged.</span>

<span class="sd">   Translated from MPFIT (Craig Markwardt&#39;s IDL package) to Python,</span>
<span class="sd">   August, 2002.  Mark Rivers</span>
<span class="sd">   Converted from Numeric to numpy (Sergey Koposov, July 2008)</span>
<span class="sd">&quot;&quot;&quot;</span>

<span class="kn">import</span> <span class="nn">numpy</span>
<span class="kn">import</span> <span class="nn">types</span>


<span class="c">#    Original FORTRAN documentation</span>
<span class="c">#    **********</span>
<span class="c">#</span>
<span class="c">#    subroutine lmdif</span>
<span class="c">#</span>
<span class="c">#    the purpose of lmdif is to minimize the sum of the squares of</span>
<span class="c">#    m nonlinear functions in n variables by a modification of</span>
<span class="c">#    the levenberg-marquardt algorithm. the user must provide a</span>
<span class="c">#    subroutine which calculates the functions. the jacobian is</span>
<span class="c">#    then calculated by a forward-difference approximation.</span>
<span class="c">#</span>
<span class="c">#    the subroutine statement is</span>
<span class="c">#</span>
<span class="c">#      subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn,</span>
<span class="c">#                       diag,mode,factor,nprint,info,nfev,fjac,</span>
<span class="c">#                       ldfjac,ipvt,qtf,wa1,wa2,wa3,wa4)</span>
<span class="c">#</span>
<span class="c">#    where</span>
<span class="c">#</span>
<span class="c">#      fcn is the name of the user-supplied subroutine which</span>
<span class="c">#        calculates the functions. fcn must be declared</span>
<span class="c">#        in an external statement in the user calling</span>
<span class="c">#        program, and should be written as follows.</span>
<span class="c">#</span>
<span class="c">#        subroutine fcn(m,n,x,fvec,iflag)</span>
<span class="c">#        integer m,n,iflag</span>
<span class="c">#        double precision x(n),fvec(m)</span>
<span class="c">#        ----------</span>
<span class="c">#        calculate the functions at x and</span>
<span class="c">#        return this vector in fvec.</span>
<span class="c">#        ----------</span>
<span class="c">#        return</span>
<span class="c">#        end</span>
<span class="c">#</span>
<span class="c">#        the value of iflag should not be changed by fcn unless</span>
<span class="c">#        the user wants to terminate execution of lmdif.</span>
<span class="c">#        in this case set iflag to a negative integer.</span>
<span class="c">#</span>
<span class="c">#      m is a positive integer input variable set to the number</span>
<span class="c">#        of functions.</span>
<span class="c">#</span>
<span class="c">#      n is a positive integer input variable set to the number</span>
<span class="c">#        of variables. n must not exceed m.</span>
<span class="c">#</span>
<span class="c">#      x is an array of length n. on input x must contain</span>
<span class="c">#        an initial estimate of the solution vector. on output x</span>
<span class="c">#        contains the final estimate of the solution vector.</span>
<span class="c">#</span>
<span class="c">#      fvec is an output array of length m which contains</span>
<span class="c">#        the functions evaluated at the output x.</span>
<span class="c">#</span>
<span class="c">#      ftol is a nonnegative input variable. termination</span>
<span class="c">#        occurs when both the actual and predicted relative</span>
<span class="c">#        reductions in the sum of squares are at most ftol.</span>
<span class="c">#        therefore, ftol measures the relative error desired</span>
<span class="c">#        in the sum of squares.</span>
<span class="c">#</span>
<span class="c">#      xtol is a nonnegative input variable. termination</span>
<span class="c">#        occurs when the relative error between two consecutive</span>
<span class="c">#        iterates is at most xtol. therefore, xtol measures the</span>
<span class="c">#        relative error desired in the approximate solution.</span>
<span class="c">#</span>
<span class="c">#      gtol is a nonnegative input variable. termination</span>
<span class="c">#        occurs when the cosine of the angle between fvec and</span>
<span class="c">#        any column of the jacobian is at most gtol in absolute</span>
<span class="c">#        value. therefore, gtol measures the orthogonality</span>
<span class="c">#        desired between the function vector and the columns</span>
<span class="c">#        of the jacobian.</span>
<span class="c">#</span>
<span class="c">#      maxfev is a positive integer input variable. termination</span>
<span class="c">#        occurs when the number of calls to fcn is at least</span>
<span class="c">#        maxfev by the end of an iteration.</span>
<span class="c">#</span>
<span class="c">#      epsfcn is an input variable used in determining a suitable</span>
<span class="c">#        step length for the forward-difference approximation. this</span>
<span class="c">#        approximation assumes that the relative errors in the</span>
<span class="c">#        functions are of the order of epsfcn. if epsfcn is less</span>
<span class="c">#        than the machine precision, it is assumed that the relative</span>
<span class="c">#        errors in the functions are of the order of the machine</span>
<span class="c">#        precision.</span>
<span class="c">#</span>
<span class="c">#      diag is an array of length n. if mode = 1 (see</span>
<span class="c">#        below), diag is internally set. if mode = 2, diag</span>
<span class="c">#        must contain positive entries that serve as</span>
<span class="c">#        multiplicative scale factors for the variables.</span>
<span class="c">#</span>
<span class="c">#      mode is an integer input variable. if mode = 1, the</span>
<span class="c">#        variables will be scaled internally. if mode = 2,</span>
<span class="c">#        the scaling is specified by the input diag. other</span>
<span class="c">#        values of mode are equivalent to mode = 1.</span>
<span class="c">#</span>
<span class="c">#      factor is a positive input variable used in determining the</span>
<span class="c">#        initial step bound. this bound is set to the product of</span>
<span class="c">#        factor and the euclidean norm of diag*x if nonzero, or else</span>
<span class="c">#        to factor itself. in most cases factor should lie in the</span>
<span class="c">#        interval (.1,100.). 100. is a generally recommended value.</span>
<span class="c">#</span>
<span class="c">#      nprint is an integer input variable that enables controlled</span>
<span class="c">#        printing of iterates if it is positive. in this case,</span>
<span class="c">#        fcn is called with iflag = 0 at the beginning of the first</span>
<span class="c">#        iteration and every nprint iterations thereafter and</span>
<span class="c">#        immediately prior to return, with x and fvec available</span>
<span class="c">#        for printing. if nprint is not positive, no special calls</span>
<span class="c">#        of fcn with iflag = 0 are made.</span>
<span class="c">#</span>
<span class="c">#      info is an integer output variable. if the user has</span>
<span class="c">#        terminated execution, info is set to the (negative)</span>
<span class="c">#        value of iflag. see description of fcn. otherwise,</span>
<span class="c">#        info is set as follows.</span>
<span class="c">#</span>
<span class="c">#        info = 0  improper input parameters.</span>
<span class="c">#</span>
<span class="c">#        info = 1  both actual and predicted relative reductions</span>
<span class="c">#                  in the sum of squares are at most ftol.</span>
<span class="c">#</span>
<span class="c">#        info = 2  relative error between two consecutive iterates</span>
<span class="c">#                  is at most xtol.</span>
<span class="c">#</span>
<span class="c">#        info = 3  conditions for info = 1 and info = 2 both hold.</span>
<span class="c">#</span>
<span class="c">#        info = 4  the cosine of the angle between fvec and any</span>
<span class="c">#                  column of the jacobian is at most gtol in</span>
<span class="c">#                  absolute value.</span>
<span class="c">#</span>
<span class="c">#        info = 5  number of calls to fcn has reached or</span>
<span class="c">#                  exceeded maxfev.</span>
<span class="c">#</span>
<span class="c">#        info = 6  ftol is too small. no further reduction in</span>
<span class="c">#                  the sum of squares is possible.</span>
<span class="c">#</span>
<span class="c">#        info = 7  xtol is too small. no further improvement in</span>
<span class="c">#                  the approximate solution x is possible.</span>
<span class="c">#</span>
<span class="c">#        info = 8  gtol is too small. fvec is orthogonal to the</span>
<span class="c">#                  columns of the jacobian to machine precision.</span>
<span class="c">#</span>
<span class="c">#      nfev is an integer output variable set to the number of</span>
<span class="c">#        calls to fcn.</span>
<span class="c">#</span>
<span class="c">#      fjac is an output m by n array. the upper n by n submatrix</span>
<span class="c">#        of fjac contains an upper triangular matrix r with</span>
<span class="c">#        diagonal elements of nonincreasing magnitude such that</span>
<span class="c">#</span>
<span class="c">#               t    t         t</span>
<span class="c">#              p *(jac *jac)*p = r *r,</span>
<span class="c">#</span>
<span class="c">#        where p is a permutation matrix and jac is the final</span>
<span class="c">#        calculated jacobian. column j of p is column ipvt(j)</span>
<span class="c">#        (see below) of the identity matrix. the lower trapezoidal</span>
<span class="c">#        part of fjac contains information generated during</span>
<span class="c">#        the computation of r.</span>
<span class="c">#</span>
<span class="c">#      ldfjac is a positive integer input variable not less than m</span>
<span class="c">#        which specifies the leading dimension of the array fjac.</span>
<span class="c">#</span>
<span class="c">#      ipvt is an integer output array of length n. ipvt</span>
<span class="c">#        defines a permutation matrix p such that jac*p = q*r,</span>
<span class="c">#        where jac is the final calculated jacobian, q is</span>
<span class="c">#        orthogonal (not stored), and r is upper triangular</span>
<span class="c">#        with diagonal elements of nonincreasing magnitude.</span>
<span class="c">#        column j of p is column ipvt(j) of the identity matrix.</span>
<span class="c">#</span>
<span class="c">#      qtf is an output array of length n which contains</span>
<span class="c">#        the first n elements of the vector (q transpose)*fvec.</span>
<span class="c">#</span>
<span class="c">#      wa1, wa2, and wa3 are work arrays of length n.</span>
<span class="c">#</span>
<span class="c">#      wa4 is a work array of length m.</span>
<span class="c">#</span>
<span class="c">#    subprograms called</span>
<span class="c">#</span>
<span class="c">#      user-supplied ...... fcn</span>
<span class="c">#</span>
<span class="c">#      minpack-supplied ... dpmpar,enorm,fdjac2,,qrfac</span>
<span class="c">#</span>
<span class="c">#      fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod</span>
<span class="c">#</span>
<span class="c">#    argonne national laboratory. minpack project. march 1980.</span>
<span class="c">#    burton s. garbow, kenneth e. hillstrom, jorge j. more</span>
<span class="c">#</span>
<span class="c">#    **********</span>

<div class="viewcode-block" id="mpfit"><a class="viewcode-back" href="../../agpy.html#agpy.mpfit_custom.mpfit">[docs]</a><span class="k">class</span> <span class="nc">mpfit</span><span class="p">:</span>
    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">fcn</span><span class="p">,</span> <span class="n">xall</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">functkw</span><span class="o">=</span><span class="p">{},</span> <span class="n">parinfo</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span>
                 <span class="n">ftol</span><span class="o">=</span><span class="mf">1.e-10</span><span class="p">,</span> <span class="n">xtol</span><span class="o">=</span><span class="mf">1.e-10</span><span class="p">,</span> <span class="n">gtol</span><span class="o">=</span><span class="mf">1.e-10</span><span class="p">,</span>
                 <span class="n">damp</span><span class="o">=</span><span class="mf">0.</span><span class="p">,</span> <span class="n">maxiter</span><span class="o">=</span><span class="mi">200</span><span class="p">,</span> <span class="n">factor</span><span class="o">=</span><span class="mf">100.</span><span class="p">,</span> <span class="n">nprint</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span>
                 <span class="n">iterfunct</span><span class="o">=</span><span class="s">&#39;default&#39;</span><span class="p">,</span> <span class="n">iterkw</span><span class="o">=</span><span class="p">{},</span> <span class="n">nocovar</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span>
                 <span class="n">fastnorm</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">rescale</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">autoderivative</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">quiet</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span>
                 <span class="n">diag</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">epsfcn</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">debug</span><span class="o">=</span><span class="mi">0</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">  Inputs:</span>
<span class="sd">    fcn:</span>
<span class="sd">       The function to be minimized.  The function should return the weighted</span>
<span class="sd">       deviations between the model and the data, as described above.</span>

<span class="sd">    xall:</span>
<span class="sd">       An array of starting values for each of the parameters of the model.</span>
<span class="sd">       The number of parameters should be fewer than the number of measurements.</span>

<span class="sd">       This parameter is optional if the parinfo keyword is used (but see</span>
<span class="sd">       parinfo).  The parinfo keyword provides a mechanism to fix or constrain</span>
<span class="sd">       individual parameters.</span>

<span class="sd">  Keywords:</span>

<span class="sd">     autoderivative:</span>
<span class="sd">        If this is set, derivatives of the function will be computed</span>
<span class="sd">        automatically via a finite differencing procedure.  If not set, then</span>
<span class="sd">        fcn must provide the (analytical) derivatives.</span>
<span class="sd">           Default: set (=1)</span>
<span class="sd">           NOTE: to supply your own analytical derivatives,</span>
<span class="sd">                 explicitly pass autoderivative=0</span>

<span class="sd">     fastnorm:</span>
<span class="sd">        Set this keyword to select a faster algorithm to compute sum-of-square</span>
<span class="sd">        values internally.  For systems with large numbers of data points, the</span>
<span class="sd">        standard algorithm can become prohibitively slow because it cannot be</span>
<span class="sd">        vectorized well.  By setting this keyword, MPFIT will run faster, but</span>
<span class="sd">        it will be more prone to floating point overflows and underflows.  Thus, setting</span>
<span class="sd">        this keyword may sacrifice some stability in the fitting process.</span>
<span class="sd">           Default: clear (=0)</span>

<span class="sd">     ftol:</span>
<span class="sd">        A nonnegative input variable. Termination occurs when both the actual</span>
<span class="sd">        and predicted relative reductions in the sum of squares are at most</span>
<span class="sd">        ftol (and status is accordingly set to 1 or 3).  Therefore, ftol</span>
<span class="sd">        measures the relative error desired in the sum of squares.</span>
<span class="sd">           Default: 1E-10</span>

<span class="sd">     functkw:</span>
<span class="sd">        A dictionary which contains the parameters to be passed to the</span>
<span class="sd">        user-supplied function specified by fcn via the standard Python</span>
<span class="sd">        keyword dictionary mechanism.  This is the way you can pass additional</span>
<span class="sd">        data to your user-supplied function without using global variables.</span>

<span class="sd">        Consider the following example:</span>
<span class="sd">           if functkw = {&#39;xval&#39;:[1.,2.,3.], &#39;yval&#39;:[1.,4.,9.],</span>
<span class="sd">                         &#39;errval&#39;:[1.,1.,1.] }</span>
<span class="sd">        then the user supplied function should be declared like this:</span>
<span class="sd">           def myfunct(p, fjac=None, xval=None, yval=None, errval=None):</span>

<span class="sd">        Default: {}   No extra parameters are passed to the user-supplied</span>
<span class="sd">                      function.</span>

<span class="sd">     gtol:</span>
<span class="sd">        A nonnegative input variable. Termination occurs when the cosine of</span>
<span class="sd">        the angle between fvec and any column of the jacobian is at most gtol</span>
<span class="sd">        in absolute value (and status is accordingly set to 4). Therefore,</span>
<span class="sd">        gtol measures the orthogonality desired between the function vector</span>
<span class="sd">        and the columns of the jacobian.</span>
<span class="sd">           Default: 1e-10</span>

<span class="sd">     iterkw:</span>
<span class="sd">        The keyword arguments to be passed to iterfunct via the dictionary</span>
<span class="sd">        keyword mechanism.  This should be a dictionary and is similar in</span>
<span class="sd">        operation to FUNCTKW.</span>
<span class="sd">           Default: {}  No arguments are passed.</span>

<span class="sd">     iterfunct:</span>
<span class="sd">        The name of a function to be called upon each NPRINT iteration of the</span>
<span class="sd">        MPFIT routine.  It should be declared in the following way:</span>
<span class="sd">           def iterfunct(myfunct, p, iter, fnorm, functkw=None,</span>
<span class="sd">                         parinfo=None, quiet=0, dof=None, [iterkw keywords here])</span>
<span class="sd">           # perform custom iteration update</span>

<span class="sd">        iterfunct must accept all three keyword parameters (FUNCTKW, PARINFO</span>
<span class="sd">        and QUIET).</span>

<span class="sd">        myfunct:  The user-supplied function to be minimized,</span>
<span class="sd">        p:      The current set of model parameters</span>
<span class="sd">        iter:    The iteration number</span>
<span class="sd">        functkw:  The arguments to be passed to myfunct.</span>
<span class="sd">        fnorm:  The chi-squared value.</span>
<span class="sd">        quiet:  Set when no textual output should be printed.</span>
<span class="sd">        dof:      The number of degrees of freedom, normally the number of points</span>
<span class="sd">                  less the number of free parameters.</span>
<span class="sd">        See below for documentation of parinfo.</span>

<span class="sd">        In implementation, iterfunct can perform updates to the terminal or</span>
<span class="sd">        graphical user interface, to provide feedback while the fit proceeds.</span>
<span class="sd">        If the fit is to be stopped for any reason, then iterfunct should return a</span>
<span class="sd">        a status value between -15 and -1.  Otherwise it should return None</span>
<span class="sd">        (e.g. no return statement) or 0.</span>
<span class="sd">        In principle, iterfunct should probably not modify the parameter values,</span>
<span class="sd">        because it may interfere with the algorithm&#39;s stability.  In practice it</span>
<span class="sd">        is allowed.</span>

<span class="sd">        Default: an internal routine is used to print the parameter values.</span>

<span class="sd">        Set iterfunct=None if there is no user-defined routine and you don&#39;t</span>
<span class="sd">        want the internal default routine be called.</span>

<span class="sd">     maxiter:</span>
<span class="sd">        The maximum number of iterations to perform.  If the number is exceeded,</span>
<span class="sd">        then the status value is set to 5 and MPFIT returns.</span>
<span class="sd">        Default: 200 iterations</span>

<span class="sd">     nocovar:</span>
<span class="sd">        Set this keyword to prevent the calculation of the covariance matrix</span>
<span class="sd">        before returning (see COVAR)</span>
<span class="sd">        Default: clear (=0)  The covariance matrix is returned</span>

<span class="sd">     nprint:</span>
<span class="sd">        The frequency with which iterfunct is called.  A value of 1 indicates</span>
<span class="sd">        that iterfunct is called with every iteration, while 2 indicates every</span>
<span class="sd">        other iteration, etc.  Note that several Levenberg-Marquardt attempts</span>
<span class="sd">        can be made in a single iteration.</span>
<span class="sd">        Default value: 1</span>

<span class="sd">     parinfo</span>
<span class="sd">        Provides a mechanism for more sophisticated constraints to be placed on</span>
<span class="sd">        parameter values.  When parinfo is not passed, then it is assumed that</span>
<span class="sd">        all parameters are free and unconstrained.  Values in parinfo are never</span>
<span class="sd">        modified during a call to MPFIT.</span>

<span class="sd">        See description above for the structure of PARINFO.</span>

<span class="sd">        Default value: None  All parameters are free and unconstrained.</span>

<span class="sd">     quiet:</span>
<span class="sd">        Set this keyword when no textual output should be printed by MPFIT</span>

<span class="sd">     damp:</span>
<span class="sd">        A scalar number, indicating the cut-off value of residuals where</span>
<span class="sd">        &quot;damping&quot; will occur.  Residuals with magnitudes greater than this</span>
<span class="sd">        number will be replaced by their hyperbolic tangent.  This partially</span>
<span class="sd">        mitigates the so-called large residual problem inherent in</span>
<span class="sd">        least-squares solvers (as for the test problem CURVI,</span>
<span class="sd">        http://www.maxthis.com/curviex.htm).</span>
<span class="sd">        A value of 0 indicates no damping.</span>
<span class="sd">           Default: 0</span>

<span class="sd">        Note: DAMP doesn&#39;t work with autoderivative=0</span>

<span class="sd">     xtol:</span>
<span class="sd">        A nonnegative input variable. Termination occurs when the relative error</span>
<span class="sd">        between two consecutive iterates is at most xtol (and status is</span>
<span class="sd">        accordingly set to 2 or 3).  Therefore, xtol measures the relative error</span>
<span class="sd">        desired in the approximate solution.</span>
<span class="sd">        Default: 1E-10</span>

<span class="sd">   Outputs:</span>

<span class="sd">     Returns an object of type mpfit.  The results are attributes of this class,</span>
<span class="sd">     e.g. mpfit.status, mpfit.errmsg, mpfit.params, npfit.niter, mpfit.covar.</span>

<span class="sd">     .status</span>
<span class="sd">        An integer status code is returned.  All values greater than zero can</span>
<span class="sd">        represent success (however .status == 5 may indicate failure to</span>
<span class="sd">        converge). It can have one of the following values:</span>

<span class="sd">        -16</span>
<span class="sd">           A parameter or function value has become infinite or an undefined</span>
<span class="sd">           number.  This is usually a consequence of numerical overflow in the</span>
<span class="sd">           user&#39;s model function, which must be avoided.</span>

<span class="sd">        -15 to -1</span>
<span class="sd">           These are error codes that either MYFUNCT or iterfunct may return to</span>
<span class="sd">           terminate the fitting process.  Values from -15 to -1 are reserved</span>
<span class="sd">           for the user functions and will not clash with MPFIT.</span>

<span class="sd">        0  Improper input parameters.</span>

<span class="sd">        1  Both actual and predicted relative reductions in the sum of squares</span>
<span class="sd">           are at most ftol.</span>

<span class="sd">        2  Relative error between two consecutive iterates is at most xtol</span>

<span class="sd">        3  Conditions for status = 1 and status = 2 both hold.</span>

<span class="sd">        4  The cosine of the angle between fvec and any column of the jacobian</span>
<span class="sd">           is at most gtol in absolute value.</span>

<span class="sd">        5  The maximum number of iterations has been reached.</span>

<span class="sd">        6  ftol is too small. No further reduction in the sum of squares is</span>
<span class="sd">           possible.</span>

<span class="sd">        7  xtol is too small. No further improvement in the approximate solution</span>
<span class="sd">           x is possible.</span>

<span class="sd">        8  gtol is too small. fvec is orthogonal to the columns of the jacobian</span>
<span class="sd">           to machine precision.</span>

<span class="sd">     .fnorm</span>
<span class="sd">        The value of the summed squared residuals for the returned parameter</span>
<span class="sd">        values.</span>

<span class="sd">     .covar</span>
<span class="sd">        The covariance matrix for the set of parameters returned by MPFIT.</span>
<span class="sd">        The matrix is NxN where N is the number of  parameters.  The square root</span>
<span class="sd">        of the diagonal elements gives the formal 1-sigma statistical errors on</span>
<span class="sd">        the parameters if errors were treated &quot;properly&quot; in fcn.</span>
<span class="sd">        Parameter errors are also returned in .perror.</span>

<span class="sd">        To compute the correlation matrix, pcor, use this example:</span>
<span class="sd">           cov = mpfit.covar</span>
<span class="sd">           pcor = cov * 0.</span>
<span class="sd">           for i in range(n):</span>
<span class="sd">              for j in range(n):</span>
<span class="sd">                 pcor[i,j] = cov[i,j]/sqrt(cov[i,i]*cov[j,j])</span>

<span class="sd">        If nocovar is set or MPFIT terminated abnormally, then .covar is set to</span>
<span class="sd">        a scalar with value None.</span>

<span class="sd">     .errmsg</span>
<span class="sd">        A string error or warning message is returned.</span>

<span class="sd">     .nfev</span>
<span class="sd">        The number of calls to MYFUNCT performed.</span>

<span class="sd">     .niter</span>
<span class="sd">        The number of iterations completed.</span>

<span class="sd">     .perror</span>
<span class="sd">        The formal 1-sigma errors in each parameter, computed from the</span>
<span class="sd">        covariance matrix.  If a parameter is held fixed, or if it touches a</span>
<span class="sd">        boundary, then the error is reported as zero.</span>

<span class="sd">        If the fit is unweighted (i.e. no errors were given, or the weights</span>
<span class="sd">        were uniformly set to unity), then .perror will probably not represent</span>
<span class="sd">        the true parameter uncertainties.</span>

<span class="sd">        *If* you can assume that the true reduced chi-squared value is unity --</span>
<span class="sd">        meaning that the fit is implicitly assumed to be of good quality --</span>
<span class="sd">        then the estimated parameter uncertainties can be computed by scaling</span>
<span class="sd">        .perror by the measured chi-squared value.</span>

<span class="sd">           dof = len(x) - len(mpfit.params) # deg of freedom</span>
<span class="sd">           # scaled uncertainties</span>
<span class="sd">           pcerror = mpfit.perror * sqrt(mpfit.fnorm / dof)</span>

<span class="sd">        &quot;&quot;&quot;</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">niter</span> <span class="o">=</span> <span class="mi">0</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">params</span> <span class="o">=</span> <span class="bp">None</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">covar</span> <span class="o">=</span> <span class="bp">None</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">perror</span> <span class="o">=</span> <span class="bp">None</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">=</span> <span class="mi">0</span>  <span class="c"># Invalid input flag set while we check inputs</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">debug</span> <span class="o">=</span> <span class="n">debug</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;&#39;</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">fastnorm</span> <span class="o">=</span> <span class="n">fastnorm</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">nfev</span> <span class="o">=</span> <span class="mi">0</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">damp</span> <span class="o">=</span> <span class="n">damp</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">machar</span> <span class="o">=</span> <span class="n">machar</span><span class="p">(</span><span class="n">double</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
        <span class="n">machep</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">machar</span><span class="o">.</span><span class="n">machep</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">dof</span><span class="o">=</span><span class="mi">0</span>

        <span class="k">if</span> <span class="p">(</span><span class="n">fcn</span><span class="o">==</span><span class="bp">None</span><span class="p">):</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&quot;Usage: parms = mpfit(&#39;myfunt&#39;, ... )&quot;</span>
            <span class="k">return</span>

        <span class="k">if</span> <span class="p">(</span><span class="n">iterfunct</span> <span class="o">==</span> <span class="s">&#39;default&#39;</span><span class="p">):</span> <span class="n">iterfunct</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">defiter</span>

        <span class="c">## Parameter damping doesn&#39;t work when user is providing their own</span>
        <span class="c">## gradients.</span>
        <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">damp</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">autoderivative</span> <span class="o">==</span> <span class="mi">0</span><span class="p">):</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span>  <span class="s">&#39;ERROR: keywords DAMP and AUTODERIVATIVE are mutually exclusive&#39;</span>
            <span class="k">return</span>

        <span class="c">## Parameters can either be stored in parinfo, or x. x takes precedence if it exists</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">xall</span> <span class="o">==</span> <span class="bp">None</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">parinfo</span> <span class="o">==</span> <span class="bp">None</span><span class="p">):</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;ERROR: must pass parameters in P or PARINFO&#39;</span>
            <span class="k">return</span>

        <span class="c">## Be sure that PARINFO is of the right type</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">parinfo</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">):</span>
            <span class="k">if</span> <span class="p">(</span><span class="nb">type</span><span class="p">(</span><span class="n">parinfo</span><span class="p">)</span> <span class="o">!=</span> <span class="n">types</span><span class="o">.</span><span class="n">ListType</span><span class="p">):</span>
                <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;ERROR: PARINFO must be a list of dictionaries.&#39;</span>
                <span class="k">return</span>
            <span class="k">else</span><span class="p">:</span>
                <span class="k">if</span> <span class="p">(</span><span class="nb">type</span><span class="p">(</span><span class="n">parinfo</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span> <span class="o">!=</span> <span class="n">types</span><span class="o">.</span><span class="n">DictionaryType</span><span class="p">):</span>
                    <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;ERROR: PARINFO must be a list of dictionaries.&#39;</span>
                    <span class="k">return</span>
            <span class="k">if</span> <span class="p">((</span><span class="n">xall</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">xall</span><span class="p">)</span> <span class="o">!=</span> <span class="nb">len</span><span class="p">(</span><span class="n">parinfo</span><span class="p">))):</span>
                <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;ERROR: number of elements in PARINFO and P must agree&#39;</span>
                <span class="k">return</span>

        <span class="c">## If the parameters were not specified at the command line, then</span>
        <span class="c">## extract them from PARINFO</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">xall</span> <span class="o">==</span> <span class="bp">None</span><span class="p">):</span>
            <span class="n">xall</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">parinfo</span><span class="p">(</span><span class="n">parinfo</span><span class="p">,</span> <span class="s">&#39;value&#39;</span><span class="p">)</span>
            <span class="k">if</span> <span class="p">(</span><span class="n">xall</span> <span class="o">==</span> <span class="bp">None</span><span class="p">):</span>
                <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;ERROR: either P or PARINFO(*)[&quot;value&quot;] must be supplied.&#39;</span>
                <span class="k">return</span>

        <span class="c">## Make sure parameters are Numeric arrays of type Float</span>
        <span class="n">xall</span> <span class="o">=</span> <span class="n">asarray</span><span class="p">(</span><span class="n">xall</span><span class="p">,</span> <span class="nb">float</span><span class="p">)</span>

        <span class="n">npar</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">xall</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span>  <span class="o">=</span> <span class="o">-</span><span class="mf">1.</span>
        <span class="n">fnorm1</span> <span class="o">=</span> <span class="o">-</span><span class="mf">1.</span>

        <span class="c">## TIED parameters?</span>
        <span class="n">ptied</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">parinfo</span><span class="p">(</span><span class="n">parinfo</span><span class="p">,</span> <span class="s">&#39;tied&#39;</span><span class="p">,</span> <span class="n">default</span><span class="o">=</span><span class="s">&#39;&#39;</span><span class="p">,</span> <span class="n">n</span><span class="o">=</span><span class="n">npar</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">qanytied</span> <span class="o">=</span> <span class="mi">0</span>
        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">npar</span><span class="p">):</span>
            <span class="n">ptied</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">ptied</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">strip</span><span class="p">()</span>
            <span class="k">if</span> <span class="p">(</span><span class="n">ptied</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">!=</span> <span class="s">&#39;&#39;</span><span class="p">):</span> <span class="bp">self</span><span class="o">.</span><span class="n">qanytied</span> <span class="o">=</span> <span class="mi">1</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">ptied</span> <span class="o">=</span> <span class="n">ptied</span>

        <span class="c">## FIXED parameters ?</span>
        <span class="n">pfixed</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">parinfo</span><span class="p">(</span><span class="n">parinfo</span><span class="p">,</span> <span class="s">&#39;fixed&#39;</span><span class="p">,</span> <span class="n">default</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">n</span><span class="o">=</span><span class="n">npar</span><span class="p">)</span>
        <span class="n">pfixed</span> <span class="o">=</span> <span class="p">(</span><span class="n">pfixed</span> <span class="o">==</span> <span class="mi">1</span><span class="p">)</span>
        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">npar</span><span class="p">):</span>
            <span class="n">pfixed</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">pfixed</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="ow">or</span> <span class="p">(</span><span class="n">ptied</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">!=</span> <span class="s">&#39;&#39;</span><span class="p">)</span> <span class="c">## Tied parameters are also effectively fixed</span>

        <span class="c">## Finite differencing step, absolute and relative, and sidedness of deriv.</span>
        <span class="n">step</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">parinfo</span><span class="p">(</span><span class="n">parinfo</span><span class="p">,</span> <span class="s">&#39;step&#39;</span><span class="p">,</span> <span class="n">default</span><span class="o">=</span><span class="mf">0.</span><span class="p">,</span> <span class="n">n</span><span class="o">=</span><span class="n">npar</span><span class="p">)</span>
        <span class="n">dstep</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">parinfo</span><span class="p">(</span><span class="n">parinfo</span><span class="p">,</span> <span class="s">&#39;relstep&#39;</span><span class="p">,</span> <span class="n">default</span><span class="o">=</span><span class="mf">0.</span><span class="p">,</span> <span class="n">n</span><span class="o">=</span><span class="n">npar</span><span class="p">)</span>
        <span class="n">dside</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">parinfo</span><span class="p">(</span><span class="n">parinfo</span><span class="p">,</span> <span class="s">&#39;mpside&#39;</span><span class="p">,</span>  <span class="n">default</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">n</span><span class="o">=</span><span class="n">npar</span><span class="p">)</span>

        <span class="c">## Maximum and minimum steps allowed to be taken in one iteration</span>
        <span class="n">maxstep</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">parinfo</span><span class="p">(</span><span class="n">parinfo</span><span class="p">,</span> <span class="s">&#39;mpmaxstep&#39;</span><span class="p">,</span> <span class="n">default</span><span class="o">=</span><span class="mf">0.</span><span class="p">,</span> <span class="n">n</span><span class="o">=</span><span class="n">npar</span><span class="p">)</span>
        <span class="n">minstep</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">parinfo</span><span class="p">(</span><span class="n">parinfo</span><span class="p">,</span> <span class="s">&#39;mpminstep&#39;</span><span class="p">,</span> <span class="n">default</span><span class="o">=</span><span class="mf">0.</span><span class="p">,</span> <span class="n">n</span><span class="o">=</span><span class="n">npar</span><span class="p">)</span>
        <span class="n">qmin</span> <span class="o">=</span> <span class="n">minstep</span> <span class="o">*</span> <span class="mi">0</span>  <span class="c">## Remove minstep for now!!</span>
        <span class="n">qmax</span> <span class="o">=</span> <span class="n">maxstep</span> <span class="o">!=</span> <span class="mi">0</span>
        <span class="n">wh</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">(((</span><span class="n">qmin</span><span class="o">!=</span><span class="mf">0.</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">qmax</span><span class="o">!=</span><span class="mf">0.</span><span class="p">))</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">maxstep</span> <span class="o">&lt;</span> <span class="n">minstep</span><span class="p">)))[</span><span class="mi">0</span><span class="p">]</span>
        <span class="k">if</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">wh</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;ERROR: MPMINSTEP is greater than MPMAXSTEP&#39;</span>
            <span class="k">return</span>
        <span class="n">wh</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">((</span><span class="n">qmin</span><span class="o">!=</span><span class="mf">0.</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">qmax</span><span class="o">!=</span><span class="mf">0.</span><span class="p">)))[</span><span class="mi">0</span><span class="p">]</span>
        <span class="n">qminmax</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">wh</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">)</span>

        <span class="c">## Finish up the free parameters</span>
        <span class="n">ifree</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">(</span><span class="n">pfixed</span> <span class="o">!=</span> <span class="mi">1</span><span class="p">))[</span><span class="mi">0</span><span class="p">]</span>
        <span class="n">nfree</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">ifree</span><span class="p">)</span>
        <span class="k">if</span> <span class="n">nfree</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;ERROR: no free parameters&#39;</span>
            <span class="k">return</span>

        <span class="c">## Compose only VARYING parameters</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">params</span> <span class="o">=</span> <span class="n">xall</span>    <span class="c">## self.params is the set of parameters to be returned</span>
        <span class="n">x</span> <span class="o">=</span> <span class="n">take</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="p">,</span> <span class="n">ifree</span><span class="p">)</span>  <span class="c">## x is the set of free parameters</span>

        <span class="c">## LIMITED parameters ?</span>
        <span class="n">limited</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">parinfo</span><span class="p">(</span><span class="n">parinfo</span><span class="p">,</span> <span class="s">&#39;limited&#39;</span><span class="p">,</span> <span class="n">default</span><span class="o">=</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">],</span> <span class="n">n</span><span class="o">=</span><span class="n">npar</span><span class="p">)</span>
        <span class="n">limits</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">parinfo</span><span class="p">(</span><span class="n">parinfo</span><span class="p">,</span> <span class="s">&#39;limits&#39;</span><span class="p">,</span> <span class="n">default</span><span class="o">=</span><span class="p">[</span><span class="mf">0.</span><span class="p">,</span><span class="mf">0.</span><span class="p">],</span> <span class="n">n</span><span class="o">=</span><span class="n">npar</span><span class="p">)</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">limited</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">limits</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">):</span>
            <span class="c">## Error checking on limits in parinfo</span>
            <span class="n">wh</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">((</span><span class="n">limited</span><span class="p">[:,</span><span class="mi">0</span><span class="p">]</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">xall</span> <span class="o">&lt;</span> <span class="n">limits</span><span class="p">[:,</span><span class="mi">0</span><span class="p">]))</span> <span class="o">|</span>
                                 <span class="p">(</span><span class="n">limited</span><span class="p">[:,</span><span class="mi">1</span><span class="p">]</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">xall</span> <span class="o">&gt;</span> <span class="n">limits</span><span class="p">[:,</span><span class="mi">1</span><span class="p">]))))[</span><span class="mi">0</span><span class="p">]</span>
            <span class="k">if</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">wh</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span>
                <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;ERROR: parameters are not within PARINFO limits&#39;</span>
                <span class="k">return</span>
            <span class="n">wh</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">((</span><span class="n">limited</span><span class="p">[:,</span><span class="mi">0</span><span class="p">]</span> <span class="o">&amp;</span> <span class="n">limited</span><span class="p">[:,</span><span class="mi">1</span><span class="p">])</span> <span class="o">&amp;</span>
                                 <span class="p">(</span><span class="n">limits</span><span class="p">[:,</span><span class="mi">0</span><span class="p">]</span> <span class="o">&gt;=</span> <span class="n">limits</span><span class="p">[:,</span><span class="mi">1</span><span class="p">])</span> <span class="o">&amp;</span>
                                 <span class="p">(</span><span class="n">pfixed</span> <span class="o">==</span> <span class="mi">0</span><span class="p">)))[</span><span class="mi">0</span><span class="p">]</span>
            <span class="k">if</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">wh</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span>
                <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;ERROR: PARINFO parameter limits are not consistent&#39;</span>
                <span class="k">return</span>

            <span class="c">## Transfer structure values to local variables</span>
            <span class="n">qulim</span> <span class="o">=</span> <span class="n">take</span><span class="p">(</span><span class="n">limited</span><span class="p">[:,</span><span class="mi">1</span><span class="p">],</span> <span class="n">ifree</span><span class="p">)</span>
            <span class="n">ulim</span>  <span class="o">=</span> <span class="n">take</span><span class="p">(</span><span class="n">limits</span> <span class="p">[:,</span><span class="mi">1</span><span class="p">],</span> <span class="n">ifree</span><span class="p">)</span>
            <span class="n">qllim</span> <span class="o">=</span> <span class="n">take</span><span class="p">(</span><span class="n">limited</span><span class="p">[:,</span><span class="mi">0</span><span class="p">],</span> <span class="n">ifree</span><span class="p">)</span>
            <span class="n">llim</span>  <span class="o">=</span> <span class="n">take</span><span class="p">(</span><span class="n">limits</span> <span class="p">[:,</span><span class="mi">0</span><span class="p">],</span> <span class="n">ifree</span><span class="p">)</span>

            <span class="n">wh</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">((</span><span class="n">qulim</span><span class="o">!=</span><span class="mf">0.</span><span class="p">)</span> <span class="o">|</span> <span class="p">(</span><span class="n">qllim</span><span class="o">!=</span><span class="mf">0.</span><span class="p">)))[</span><span class="mi">0</span><span class="p">]</span>
            <span class="k">if</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">wh</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span> <span class="n">qanylim</span> <span class="o">=</span> <span class="mi">1</span>
            <span class="k">else</span><span class="p">:</span> <span class="n">qanylim</span> <span class="o">=</span> <span class="mi">0</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="c">## Fill in local variables with dummy values</span>
            <span class="n">qulim</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">nfree</span><span class="p">)</span>
            <span class="n">ulim</span>  <span class="o">=</span> <span class="n">x</span> <span class="o">*</span> <span class="mf">0.</span>
            <span class="n">qllim</span> <span class="o">=</span> <span class="n">qulim</span>
            <span class="n">llim</span>  <span class="o">=</span> <span class="n">x</span> <span class="o">*</span> <span class="mf">0.</span>
            <span class="n">qanylim</span> <span class="o">=</span> <span class="mi">0</span>

        <span class="n">n</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
        <span class="c">## Check input parameters for errors</span>
        <span class="k">if</span> <span class="p">((</span><span class="n">n</span> <span class="o">&lt;</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">or</span> <span class="p">(</span><span class="n">ftol</span> <span class="o">&lt;=</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">or</span> <span class="p">(</span><span class="n">xtol</span> <span class="o">&lt;=</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">or</span> <span class="p">(</span><span class="n">gtol</span> <span class="o">&lt;=</span> <span class="mi">0</span><span class="p">)</span>
                    <span class="ow">or</span> <span class="p">(</span><span class="n">maxiter</span> <span class="o">&lt;=</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">or</span> <span class="p">(</span><span class="n">factor</span> <span class="o">&lt;=</span> <span class="mi">0</span><span class="p">)):</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;ERROR: input keywords are inconsistent&#39;</span>
            <span class="k">return</span>

        <span class="k">if</span> <span class="p">(</span><span class="n">rescale</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">):</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;ERROR: DIAG parameter scales are inconsistent&#39;</span>
            <span class="k">if</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">diag</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">n</span><span class="p">):</span> <span class="k">return</span>
            <span class="n">wh</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">(</span><span class="n">diag</span> <span class="o">&lt;=</span> <span class="mi">0</span><span class="p">))[</span><span class="mi">0</span><span class="p">]</span>
            <span class="k">if</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">wh</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span> <span class="k">return</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;&#39;</span>

        <span class="c"># Make sure x is a Numeric array of type Float</span>
        <span class="n">x</span> <span class="o">=</span> <span class="n">asarray</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="nb">float</span><span class="p">)</span>

        <span class="p">[</span><span class="bp">self</span><span class="o">.</span><span class="n">status</span><span class="p">,</span> <span class="n">fvec</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">call</span><span class="p">(</span><span class="n">fcn</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="p">,</span> <span class="n">functkw</span><span class="p">)</span>
        <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">&lt;</span> <span class="mi">0</span><span class="p">):</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;ERROR: first call to &quot;&#39;</span><span class="o">+</span><span class="nb">str</span><span class="p">(</span><span class="n">fcn</span><span class="p">)</span><span class="o">+</span><span class="s">&#39;&quot; failed&#39;</span>
            <span class="k">return</span>

        <span class="n">m</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">fvec</span><span class="p">)</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">m</span> <span class="o">&lt;</span> <span class="n">n</span><span class="p">):</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;ERROR: number of parameters must not exceed data&#39;</span>
            <span class="k">return</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">dof</span> <span class="o">=</span> <span class="n">m</span><span class="o">-</span><span class="n">nfree</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">enorm</span><span class="p">(</span><span class="n">fvec</span><span class="p">)</span>

        <span class="c">## Initialize Levelberg-Marquardt parameter and iteration counter</span>

        <span class="n">par</span> <span class="o">=</span> <span class="mf">0.</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">niter</span> <span class="o">=</span> <span class="mi">1</span>
        <span class="n">qtf</span> <span class="o">=</span> <span class="n">x</span> <span class="o">*</span> <span class="mf">0.</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">=</span> <span class="mi">0</span>

        <span class="c">## Beginning of the outer loop</span>

        <span class="k">while</span><span class="p">(</span><span class="mi">1</span><span class="p">):</span>

            <span class="c">## If requested, call fcn to enable printing of iterates</span>
            <span class="n">put</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="p">,</span> <span class="n">ifree</span><span class="p">,</span> <span class="n">x</span><span class="p">)</span>
            <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">qanytied</span><span class="p">):</span> <span class="bp">self</span><span class="o">.</span><span class="n">params</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">tie</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="p">,</span> <span class="n">ptied</span><span class="p">)</span>

            <span class="k">if</span> <span class="p">(</span><span class="n">nprint</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">iterfunct</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">):</span>
                <span class="k">if</span> <span class="p">(((</span><span class="bp">self</span><span class="o">.</span><span class="n">niter</span><span class="o">-</span><span class="mi">1</span><span class="p">)</span> <span class="o">%</span> <span class="n">nprint</span><span class="p">)</span> <span class="o">==</span> <span class="mi">0</span><span class="p">):</span>
                    <span class="n">mperr</span> <span class="o">=</span> <span class="mi">0</span>
                    <span class="n">xnew0</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>

                    <span class="n">dof</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">max</span><span class="p">([</span><span class="nb">len</span><span class="p">(</span><span class="n">fvec</span><span class="p">)</span> <span class="o">-</span> <span class="nb">len</span><span class="p">(</span><span class="n">x</span><span class="p">),</span> <span class="mi">0</span><span class="p">])</span>
                    <span class="n">status</span> <span class="o">=</span> <span class="n">iterfunct</span><span class="p">(</span><span class="n">fcn</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">niter</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span><span class="o">**</span><span class="mi">2</span><span class="p">,</span>
                       <span class="n">functkw</span><span class="o">=</span><span class="n">functkw</span><span class="p">,</span> <span class="n">parinfo</span><span class="o">=</span><span class="n">parinfo</span><span class="p">,</span> <span class="n">quiet</span><span class="o">=</span><span class="n">quiet</span><span class="p">,</span>
                       <span class="n">dof</span><span class="o">=</span><span class="n">dof</span><span class="p">,</span> <span class="o">**</span><span class="n">iterkw</span><span class="p">)</span>
                    <span class="k">if</span> <span class="p">(</span><span class="n">status</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">):</span> <span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">=</span> <span class="n">status</span>

                    <span class="c">## Check for user termination</span>
                    <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">&lt;</span> <span class="mi">0</span><span class="p">):</span>
                        <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;WARNING: premature termination by &#39;</span> <span class="o">+</span> <span class="nb">str</span><span class="p">(</span><span class="n">iterfunct</span><span class="p">)</span>
                        <span class="k">return</span>

                    <span class="c">## If parameters were changed (grrr..) then re-tie</span>
                    <span class="k">if</span> <span class="p">(</span><span class="n">numpy</span><span class="o">.</span><span class="n">max</span><span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">xnew0</span><span class="o">-</span><span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="p">))</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span>
                        <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">qanytied</span><span class="p">):</span> <span class="bp">self</span><span class="o">.</span><span class="n">params</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">tie</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="p">,</span> <span class="n">ptied</span><span class="p">)</span>
                        <span class="n">x</span> <span class="o">=</span> <span class="n">take</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="p">,</span> <span class="n">ifree</span><span class="p">)</span>


            <span class="c">## Calculate the jacobian matrix</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">=</span> <span class="mi">2</span>
            <span class="n">catch_msg</span> <span class="o">=</span> <span class="s">&#39;calling MPFIT_FDJAC2&#39;</span>
            <span class="n">fjac</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">fdjac2</span><span class="p">(</span><span class="n">fcn</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">fvec</span><span class="p">,</span> <span class="n">step</span><span class="p">,</span> <span class="n">qulim</span><span class="p">,</span> <span class="n">ulim</span><span class="p">,</span> <span class="n">dside</span><span class="p">,</span>
                          <span class="n">epsfcn</span><span class="o">=</span><span class="n">epsfcn</span><span class="p">,</span>
                          <span class="n">autoderivative</span><span class="o">=</span><span class="n">autoderivative</span><span class="p">,</span> <span class="n">dstep</span><span class="o">=</span><span class="n">dstep</span><span class="p">,</span>
                          <span class="n">functkw</span><span class="o">=</span><span class="n">functkw</span><span class="p">,</span> <span class="n">ifree</span><span class="o">=</span><span class="n">ifree</span><span class="p">,</span> <span class="n">xall</span><span class="o">=</span><span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="p">)</span>
            <span class="k">if</span> <span class="p">(</span><span class="n">fjac</span> <span class="o">==</span> <span class="bp">None</span><span class="p">):</span>
                <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;WARNING: premature termination by FDJAC2&#39;</span>
                <span class="k">return</span>

            <span class="c">## Determine if any of the parameters are pegged at the limits</span>
            <span class="k">if</span> <span class="p">(</span><span class="n">qanylim</span><span class="p">):</span>
                <span class="n">catch_msg</span> <span class="o">=</span> <span class="s">&#39;zeroing derivatives of pegged parameters&#39;</span>
                <span class="n">whlpeg</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">(</span><span class="n">qllim</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">x</span> <span class="o">==</span> <span class="n">llim</span><span class="p">)))[</span><span class="mi">0</span><span class="p">]</span>
                <span class="n">nlpeg</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">whlpeg</span><span class="p">)</span>
                <span class="n">whupeg</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">(</span><span class="n">qulim</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">x</span> <span class="o">==</span> <span class="n">ulim</span><span class="p">)))[</span><span class="mi">0</span><span class="p">]</span>
                <span class="n">nupeg</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">whupeg</span><span class="p">)</span>
                <span class="c">## See if any &quot;pegged&quot; values should keep their derivatives</span>
                <span class="k">if</span> <span class="p">(</span><span class="n">nlpeg</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span>
                    <span class="c">## Total derivative of sum wrt lower pegged parameters</span>
                    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nlpeg</span><span class="p">):</span>
                        <span class="n">sum0</span> <span class="o">=</span> <span class="nb">sum</span><span class="p">(</span><span class="n">fvec</span> <span class="o">*</span> <span class="n">fjac</span><span class="p">[:,</span><span class="n">whlpeg</span><span class="p">[</span><span class="n">i</span><span class="p">]])</span>
                        <span class="k">if</span> <span class="p">(</span><span class="n">sum0</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span> <span class="n">fjac</span><span class="p">[:,</span><span class="n">whlpeg</span><span class="p">[</span><span class="n">i</span><span class="p">]]</span> <span class="o">=</span> <span class="mi">0</span>
                <span class="k">if</span> <span class="p">(</span><span class="n">nupeg</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span>
                    <span class="c">## Total derivative of sum wrt upper pegged parameters</span>
                    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nupeg</span><span class="p">):</span>
                        <span class="n">sum0</span> <span class="o">=</span> <span class="nb">sum</span><span class="p">(</span><span class="n">fvec</span> <span class="o">*</span> <span class="n">fjac</span><span class="p">[:,</span><span class="n">whupeg</span><span class="p">[</span><span class="n">i</span><span class="p">]])</span>
                        <span class="k">if</span> <span class="p">(</span><span class="n">sum0</span> <span class="o">&lt;</span> <span class="mi">0</span><span class="p">):</span> <span class="n">fjac</span><span class="p">[:,</span><span class="n">whupeg</span><span class="p">[</span><span class="n">i</span><span class="p">]]</span> <span class="o">=</span> <span class="mi">0</span>

            <span class="c">## Compute the QR factorization of the jacobian</span>
            <span class="p">[</span><span class="n">fjac</span><span class="p">,</span> <span class="n">ipvt</span><span class="p">,</span> <span class="n">wa1</span><span class="p">,</span> <span class="n">wa2</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">qrfac</span><span class="p">(</span><span class="n">fjac</span><span class="p">,</span> <span class="n">pivot</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
            
            <span class="c">## On the first iteration if &quot;diag&quot; is unspecified, scale</span>
            <span class="c">## according to the norms of the columns of the initial jacobian</span>
            <span class="n">catch_msg</span> <span class="o">=</span> <span class="s">&#39;rescaling diagonal elements&#39;</span>
            <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">niter</span> <span class="o">==</span> <span class="mi">1</span><span class="p">):</span>
                <span class="k">if</span> <span class="p">((</span><span class="n">rescale</span><span class="o">==</span><span class="mi">0</span><span class="p">)</span> <span class="ow">or</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">diag</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">n</span><span class="p">)):</span>
                    <span class="n">diag</span> <span class="o">=</span> <span class="n">wa2</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
                    <span class="n">wh</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">(</span><span class="n">diag</span> <span class="o">==</span> <span class="mi">0</span><span class="p">))[</span><span class="mi">0</span><span class="p">]</span>
                    <span class="n">put</span><span class="p">(</span><span class="n">diag</span><span class="p">,</span> <span class="n">wh</span><span class="p">,</span> <span class="mf">1.</span><span class="p">)</span>

                <span class="c">## On the first iteration, calculate the norm of the scaled x</span>
                <span class="c">## and initialize the step bound delta</span>
                <span class="n">wa3</span> <span class="o">=</span> <span class="n">diag</span> <span class="o">*</span> <span class="n">x</span>
                <span class="n">xnorm</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">enorm</span><span class="p">(</span><span class="n">wa3</span><span class="p">)</span>
                <span class="n">delta</span> <span class="o">=</span> <span class="n">factor</span><span class="o">*</span><span class="n">xnorm</span>
                <span class="k">if</span> <span class="p">(</span><span class="n">delta</span> <span class="o">==</span> <span class="mf">0.</span><span class="p">):</span> <span class="n">delta</span> <span class="o">=</span> <span class="n">factor</span>

            <span class="c">## Form (q transpose)*fvec and store the first n components in qtf</span>
            <span class="n">catch_msg</span> <span class="o">=</span> <span class="s">&#39;forming (q transpose)*fvec&#39;</span>
            <span class="n">wa4</span> <span class="o">=</span> <span class="n">fvec</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
            <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
                <span class="n">lj</span> <span class="o">=</span> <span class="n">ipvt</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
                <span class="n">temp3</span> <span class="o">=</span> <span class="n">fjac</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">lj</span><span class="p">]</span>
                <span class="k">if</span> <span class="p">(</span><span class="n">temp3</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">):</span>
                    <span class="n">fj</span> <span class="o">=</span> <span class="n">fjac</span><span class="p">[</span><span class="n">j</span><span class="p">:,</span><span class="n">lj</span><span class="p">]</span>
                    <span class="n">wj</span> <span class="o">=</span> <span class="n">wa4</span><span class="p">[</span><span class="n">j</span><span class="p">:]</span>
                    <span class="c">## *** optimization wa4(j:*)</span>
                    <span class="n">wa4</span><span class="p">[</span><span class="n">j</span><span class="p">:]</span> <span class="o">=</span> <span class="n">wj</span> <span class="o">-</span> <span class="n">fj</span> <span class="o">*</span> <span class="nb">sum</span><span class="p">(</span><span class="n">fj</span><span class="o">*</span><span class="n">wj</span><span class="p">)</span> <span class="o">/</span> <span class="n">temp3</span>
                <span class="n">fjac</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">lj</span><span class="p">]</span> <span class="o">=</span> <span class="n">wa1</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
                <span class="n">qtf</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">wa4</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
            <span class="c">## From this point on, only the square matrix, consisting of the</span>
            <span class="c">## triangle of R, is needed.</span>
            <span class="n">fjac</span> <span class="o">=</span> <span class="n">fjac</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">n</span><span class="p">,</span> <span class="mi">0</span><span class="p">:</span><span class="n">n</span><span class="p">]</span>
            <span class="n">fjac</span><span class="o">.</span><span class="n">shape</span> <span class="o">=</span> <span class="p">[</span><span class="n">n</span><span class="p">,</span> <span class="n">n</span><span class="p">]</span>
            <span class="n">temp</span> <span class="o">=</span> <span class="n">fjac</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
            <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
                <span class="n">temp</span><span class="p">[:,</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">fjac</span><span class="p">[:,</span> <span class="n">ipvt</span><span class="p">[</span><span class="n">i</span><span class="p">]]</span>
            <span class="n">fjac</span> <span class="o">=</span> <span class="n">temp</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>

            <span class="c">## Check for overflow.  This should be a cheap test here since FJAC</span>
            <span class="c">## has been reduced to a (small) square matrix, and the test is</span>
            <span class="c">## O(N^2).</span>
            <span class="c">#wh = where(finite(fjac) EQ 0, ct)</span>
            <span class="c">#if ct GT 0 then goto, FAIL_OVERFLOW</span>

            <span class="c">## Compute the norm of the scaled gradient</span>
            <span class="n">catch_msg</span> <span class="o">=</span> <span class="s">&#39;computing the scaled gradient&#39;</span>
            <span class="n">gnorm</span> <span class="o">=</span> <span class="mf">0.</span>
            <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">):</span>
                <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
                    <span class="n">l</span> <span class="o">=</span> <span class="n">ipvt</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
                    <span class="k">if</span> <span class="p">(</span><span class="n">wa2</span><span class="p">[</span><span class="n">l</span><span class="p">]</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">):</span>
                        <span class="n">sum0</span> <span class="o">=</span> <span class="nb">sum</span><span class="p">(</span><span class="n">fjac</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span><span class="n">j</span><span class="p">]</span><span class="o">*</span><span class="n">qtf</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">])</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span>
                        <span class="n">gnorm</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">max</span><span class="p">([</span><span class="n">gnorm</span><span class="p">,</span><span class="nb">abs</span><span class="p">(</span><span class="n">sum0</span><span class="o">/</span><span class="n">wa2</span><span class="p">[</span><span class="n">l</span><span class="p">])])</span>

            <span class="c">## Test for convergence of the gradient norm</span>
            <span class="k">if</span> <span class="p">(</span><span class="n">gnorm</span> <span class="o">&lt;=</span> <span class="n">gtol</span><span class="p">):</span>
                <span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">=</span> <span class="mi">4</span>
                <span class="k">break</span>
            <span class="k">if</span> <span class="n">maxiter</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
                <span class="k">break</span>

            <span class="c">## Rescale if necessary</span>
            <span class="k">if</span> <span class="p">(</span><span class="n">rescale</span> <span class="o">==</span> <span class="mi">0</span><span class="p">):</span>
                <span class="n">diag</span> <span class="o">=</span> <span class="n">choose</span><span class="p">(</span><span class="n">diag</span><span class="o">&gt;</span><span class="n">wa2</span><span class="p">,</span> <span class="p">(</span><span class="n">wa2</span><span class="p">,</span> <span class="n">diag</span><span class="p">))</span>

            <span class="c">## Beginning of the inner loop</span>
            <span class="k">while</span><span class="p">(</span><span class="mi">1</span><span class="p">):</span>

                <span class="c">## Determine the levenberg-marquardt parameter</span>
                <span class="n">catch_msg</span> <span class="o">=</span> <span class="s">&#39;calculating LM parameter (MPFIT_)&#39;</span>
                <span class="p">[</span><span class="n">fjac</span><span class="p">,</span> <span class="n">par</span><span class="p">,</span> <span class="n">wa1</span><span class="p">,</span> <span class="n">wa2</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">lmpar</span><span class="p">(</span><span class="n">fjac</span><span class="p">,</span> <span class="n">ipvt</span><span class="p">,</span> <span class="n">diag</span><span class="p">,</span> <span class="n">qtf</span><span class="p">,</span>
                                                     <span class="n">delta</span><span class="p">,</span> <span class="n">wa1</span><span class="p">,</span> <span class="n">wa2</span><span class="p">,</span> <span class="n">par</span><span class="o">=</span><span class="n">par</span><span class="p">)</span>
                <span class="c">## Store the direction p and x+p. Calculate the norm of p</span>
                <span class="n">wa1</span> <span class="o">=</span> <span class="o">-</span><span class="n">wa1</span>

                <span class="k">if</span> <span class="p">(</span><span class="n">qanylim</span> <span class="o">==</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">qminmax</span> <span class="o">==</span> <span class="mi">0</span><span class="p">):</span>
                    <span class="c">## No parameter limits, so just move to new position WA2</span>
                    <span class="n">alpha</span> <span class="o">=</span> <span class="mf">1.</span>
                    <span class="n">wa2</span> <span class="o">=</span> <span class="n">x</span> <span class="o">+</span> <span class="n">wa1</span>

                <span class="k">else</span><span class="p">:</span>

                    <span class="c">## Respect the limits.  If a step were to go out of bounds, then</span>
                    <span class="c">## we should take a step in the same direction but shorter distance.</span>
                    <span class="c">## The step should take us right to the limit in that case.</span>
                    <span class="n">alpha</span> <span class="o">=</span> <span class="mf">1.</span>

                    <span class="k">if</span> <span class="p">(</span><span class="n">qanylim</span><span class="p">):</span>
                        <span class="c">## Do not allow any steps out of bounds</span>
                        <span class="n">catch_msg</span> <span class="o">=</span> <span class="s">&#39;checking for a step out of bounds&#39;</span>
                        <span class="k">if</span> <span class="p">(</span><span class="n">nlpeg</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span>
                            <span class="n">put</span><span class="p">(</span><span class="n">wa1</span><span class="p">,</span> <span class="n">whlpeg</span><span class="p">,</span> <span class="n">clip</span><span class="p">(</span>
                               <span class="n">take</span><span class="p">(</span><span class="n">wa1</span><span class="p">,</span> <span class="n">whlpeg</span><span class="p">),</span> <span class="mf">0.</span><span class="p">,</span> <span class="n">numpy</span><span class="o">.</span><span class="n">max</span><span class="p">(</span><span class="n">wa1</span><span class="p">)))</span>
                        <span class="k">if</span> <span class="p">(</span><span class="n">nupeg</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span>
                            <span class="n">put</span><span class="p">(</span><span class="n">wa1</span><span class="p">,</span> <span class="n">whupeg</span><span class="p">,</span> <span class="n">clip</span><span class="p">(</span>
                               <span class="n">take</span><span class="p">(</span><span class="n">wa1</span><span class="p">,</span> <span class="n">whupeg</span><span class="p">),</span> <span class="n">numpy</span><span class="o">.</span><span class="n">min</span><span class="p">(</span><span class="n">wa1</span><span class="p">),</span> <span class="mf">0.</span><span class="p">))</span>

                        <span class="n">dwa1</span> <span class="o">=</span> <span class="nb">abs</span><span class="p">(</span><span class="n">wa1</span><span class="p">)</span> <span class="o">&gt;</span> <span class="n">machep</span>
                        <span class="n">whl</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">(((</span><span class="n">dwa1</span><span class="o">!=</span><span class="mf">0.</span><span class="p">)</span> <span class="o">&amp;</span> <span class="n">qllim</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">((</span><span class="n">x</span> <span class="o">+</span> <span class="n">wa1</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">llim</span><span class="p">)))[</span><span class="mi">0</span><span class="p">]</span>
                        <span class="k">if</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">whl</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span>
                            <span class="n">t</span> <span class="o">=</span> <span class="p">((</span><span class="n">take</span><span class="p">(</span><span class="n">llim</span><span class="p">,</span> <span class="n">whl</span><span class="p">)</span> <span class="o">-</span> <span class="n">take</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">whl</span><span class="p">))</span> <span class="o">/</span>
                                  <span class="n">take</span><span class="p">(</span><span class="n">wa1</span><span class="p">,</span> <span class="n">whl</span><span class="p">))</span>
                            <span class="n">alpha</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">min</span><span class="p">([</span><span class="n">alpha</span><span class="p">,</span> <span class="n">numpy</span><span class="o">.</span><span class="n">min</span><span class="p">(</span><span class="n">t</span><span class="p">)])</span>
                        <span class="n">whu</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">(((</span><span class="n">dwa1</span><span class="o">!=</span><span class="mf">0.</span><span class="p">)</span> <span class="o">&amp;</span> <span class="n">qulim</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">((</span><span class="n">x</span> <span class="o">+</span> <span class="n">wa1</span><span class="p">)</span> <span class="o">&gt;</span> <span class="n">ulim</span><span class="p">)))[</span><span class="mi">0</span><span class="p">]</span>
                        <span class="k">if</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">whu</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span>
                            <span class="n">t</span> <span class="o">=</span> <span class="p">((</span><span class="n">take</span><span class="p">(</span><span class="n">ulim</span><span class="p">,</span> <span class="n">whu</span><span class="p">)</span> <span class="o">-</span> <span class="n">take</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">whu</span><span class="p">))</span> <span class="o">/</span>
                                  <span class="n">take</span><span class="p">(</span><span class="n">wa1</span><span class="p">,</span> <span class="n">whu</span><span class="p">))</span>
                            <span class="n">alpha</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">min</span><span class="p">([</span><span class="n">alpha</span><span class="p">,</span> <span class="n">numpy</span><span class="o">.</span><span class="n">min</span><span class="p">(</span><span class="n">t</span><span class="p">)])</span>

                    <span class="c">## Obey any max step values.</span>
                    <span class="k">if</span> <span class="p">(</span><span class="n">qminmax</span><span class="p">):</span>
                        <span class="n">nwa1</span> <span class="o">=</span> <span class="n">wa1</span> <span class="o">*</span> <span class="n">alpha</span>
                        <span class="n">whmax</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">((</span><span class="n">qmax</span> <span class="o">!=</span> <span class="mf">0.</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">maxstep</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">)))[</span><span class="mi">0</span><span class="p">]</span>
                        <span class="k">if</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">whmax</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span>
                            <span class="n">mrat</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">max</span><span class="p">(</span><span class="n">take</span><span class="p">(</span><span class="n">nwa1</span><span class="p">,</span> <span class="n">whmax</span><span class="p">)</span> <span class="o">/</span>
                                       <span class="n">take</span><span class="p">(</span><span class="n">maxstep</span><span class="p">,</span> <span class="n">whmax</span><span class="p">))</span>
                            <span class="k">if</span> <span class="p">(</span><span class="n">mrat</span> <span class="o">&gt;</span> <span class="mi">1</span><span class="p">):</span> <span class="n">alpha</span> <span class="o">=</span> <span class="n">alpha</span> <span class="o">/</span> <span class="n">mrat</span>

                    <span class="c">## Scale the resulting vector</span>
                    <span class="n">wa1</span> <span class="o">=</span> <span class="n">wa1</span> <span class="o">*</span> <span class="n">alpha</span>
                    <span class="n">wa2</span> <span class="o">=</span> <span class="n">x</span> <span class="o">+</span> <span class="n">wa1</span>

                    <span class="c">## Adjust the final output values.  If the step put us exactly</span>
                    <span class="c">## on a boundary, make sure it is exact.</span>
                    <span class="n">wh</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">((</span><span class="n">qulim</span><span class="o">!=</span><span class="mf">0.</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">wa2</span> <span class="o">&gt;=</span> <span class="n">ulim</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">machep</span><span class="p">))))[</span><span class="mi">0</span><span class="p">]</span>
                    <span class="k">if</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">wh</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span> <span class="n">put</span><span class="p">(</span><span class="n">wa2</span><span class="p">,</span> <span class="n">wh</span><span class="p">,</span> <span class="n">take</span><span class="p">(</span><span class="n">ulim</span><span class="p">,</span> <span class="n">wh</span><span class="p">))</span>
                    <span class="n">wh</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">((</span><span class="n">qllim</span><span class="o">!=</span><span class="mf">0.</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">wa2</span> <span class="o">&lt;=</span> <span class="n">llim</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="n">machep</span><span class="p">))))[</span><span class="mi">0</span><span class="p">]</span>
                    <span class="k">if</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">wh</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span> <span class="n">put</span><span class="p">(</span><span class="n">wa2</span><span class="p">,</span> <span class="n">wh</span><span class="p">,</span> <span class="n">take</span><span class="p">(</span><span class="n">llim</span><span class="p">,</span> <span class="n">wh</span><span class="p">))</span>
                <span class="c"># endelse</span>
                <span class="n">wa3</span> <span class="o">=</span> <span class="n">diag</span> <span class="o">*</span> <span class="n">wa1</span>
                <span class="n">pnorm</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">enorm</span><span class="p">(</span><span class="n">wa3</span><span class="p">)</span>
                
                <span class="c">## On the first iteration, adjust the initial step bound</span>
                <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">niter</span> <span class="o">==</span> <span class="mi">1</span><span class="p">):</span> <span class="n">delta</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">min</span><span class="p">([</span><span class="n">delta</span><span class="p">,</span><span class="n">pnorm</span><span class="p">])</span>

                <span class="n">put</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="p">,</span> <span class="n">ifree</span><span class="p">,</span> <span class="n">wa2</span><span class="p">)</span>

                <span class="c">## Evaluate the function at x+p and calculate its norm</span>
                <span class="n">mperr</span> <span class="o">=</span> <span class="mi">0</span>
                <span class="n">catch_msg</span> <span class="o">=</span> <span class="s">&#39;calling &#39;</span><span class="o">+</span><span class="nb">str</span><span class="p">(</span><span class="n">fcn</span><span class="p">)</span>
                <span class="p">[</span><span class="bp">self</span><span class="o">.</span><span class="n">status</span><span class="p">,</span> <span class="n">wa4</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">call</span><span class="p">(</span><span class="n">fcn</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="p">,</span> <span class="n">functkw</span><span class="p">)</span>
                <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">&lt;</span> <span class="mi">0</span><span class="p">):</span>
                    <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span> <span class="o">=</span> <span class="s">&#39;WARNING: premature termination by &quot;&#39;</span><span class="o">+</span><span class="n">fcn</span><span class="o">+</span><span class="s">&#39;&quot;&#39;</span>
                    <span class="k">return</span>
                <span class="n">fnorm1</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">enorm</span><span class="p">(</span><span class="n">wa4</span><span class="p">)</span>

                <span class="c">## Compute the scaled actual reduction</span>
                <span class="n">catch_msg</span> <span class="o">=</span> <span class="s">&#39;computing convergence criteria&#39;</span>
                <span class="n">actred</span> <span class="o">=</span> <span class="o">-</span><span class="mf">1.</span>
                <span class="k">if</span> <span class="p">((</span><span class="mf">0.1</span> <span class="o">*</span> <span class="n">fnorm1</span><span class="p">)</span> <span class="o">&lt;</span> <span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span><span class="p">):</span> <span class="n">actred</span> <span class="o">=</span> <span class="o">-</span> <span class="p">(</span><span class="n">fnorm1</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> <span class="mf">1.</span>

                <span class="c">## Compute the scaled predicted reduction and the scaled directional</span>
                <span class="c">## derivative</span>
                <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
                    <span class="n">wa3</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="mi">0</span>
                    <span class="n">wa3</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">wa3</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">+</span> <span class="n">fjac</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span><span class="n">j</span><span class="p">]</span><span class="o">*</span><span class="n">wa1</span><span class="p">[</span><span class="n">ipvt</span><span class="p">[</span><span class="n">j</span><span class="p">]]</span>

                <span class="c">## Remember, alpha is the fraction of the full LM step actually</span>
                <span class="c">## taken</span>
                <span class="n">temp1</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">enorm</span><span class="p">(</span><span class="n">alpha</span><span class="o">*</span><span class="n">wa3</span><span class="p">)</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span>
                <span class="n">temp2</span> <span class="o">=</span> <span class="p">(</span><span class="n">sqrt</span><span class="p">(</span><span class="n">alpha</span><span class="o">*</span><span class="n">par</span><span class="p">)</span><span class="o">*</span><span class="n">pnorm</span><span class="p">)</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span>
                <span class="n">prered</span> <span class="o">=</span> <span class="n">temp1</span><span class="o">*</span><span class="n">temp1</span> <span class="o">+</span> <span class="p">(</span><span class="n">temp2</span><span class="o">*</span><span class="n">temp2</span><span class="p">)</span><span class="o">/</span><span class="mf">0.5</span>
                <span class="n">dirder</span> <span class="o">=</span> <span class="o">-</span><span class="p">(</span><span class="n">temp1</span><span class="o">*</span><span class="n">temp1</span> <span class="o">+</span> <span class="n">temp2</span><span class="o">*</span><span class="n">temp2</span><span class="p">)</span>
                
                <span class="c">## Compute the ratio of the actual to the predicted reduction.</span>
                <span class="n">ratio</span> <span class="o">=</span> <span class="mf">0.</span>
                <span class="k">if</span> <span class="p">(</span><span class="n">prered</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">):</span> <span class="n">ratio</span> <span class="o">=</span> <span class="n">actred</span><span class="o">/</span><span class="n">prered</span>

                <span class="c">## Update the step bound</span>
                <span class="k">if</span> <span class="p">(</span><span class="n">ratio</span> <span class="o">&lt;=</span> <span class="mf">0.25</span><span class="p">):</span>
                    <span class="k">if</span> <span class="p">(</span><span class="n">actred</span> <span class="o">&gt;=</span> <span class="mi">0</span><span class="p">):</span> <span class="n">temp</span> <span class="o">=</span> <span class="o">.</span><span class="mi">5</span>
                    <span class="k">else</span><span class="p">:</span> <span class="n">temp</span> <span class="o">=</span> <span class="o">.</span><span class="mi">5</span><span class="o">*</span><span class="n">dirder</span><span class="o">/</span><span class="p">(</span><span class="n">dirder</span> <span class="o">+</span> <span class="o">.</span><span class="mi">5</span><span class="o">*</span><span class="n">actred</span><span class="p">)</span>
                    <span class="k">if</span> <span class="p">((</span><span class="mf">0.1</span><span class="o">*</span><span class="n">fnorm1</span><span class="p">)</span> <span class="o">&gt;=</span> <span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span><span class="p">)</span> <span class="ow">or</span> <span class="p">(</span><span class="n">temp</span> <span class="o">&lt;</span> <span class="mf">0.1</span><span class="p">):</span> <span class="n">temp</span> <span class="o">=</span> <span class="mf">0.1</span>
                    <span class="n">delta</span> <span class="o">=</span> <span class="n">temp</span><span class="o">*</span><span class="n">numpy</span><span class="o">.</span><span class="n">min</span><span class="p">([</span><span class="n">delta</span><span class="p">,</span><span class="n">pnorm</span><span class="o">/</span><span class="mf">0.1</span><span class="p">])</span>
                    <span class="n">par</span> <span class="o">=</span> <span class="n">par</span><span class="o">/</span><span class="n">temp</span>
                <span class="k">else</span><span class="p">:</span>
                    <span class="k">if</span> <span class="p">(</span><span class="n">par</span> <span class="o">==</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">or</span> <span class="p">(</span><span class="n">ratio</span> <span class="o">&gt;=</span> <span class="mf">0.75</span><span class="p">):</span>
                        <span class="n">delta</span> <span class="o">=</span> <span class="n">pnorm</span><span class="o">/.</span><span class="mi">5</span>
                        <span class="n">par</span> <span class="o">=</span> <span class="o">.</span><span class="mi">5</span><span class="o">*</span><span class="n">par</span>

                <span class="c">## Test for successful iteration</span>
                <span class="k">if</span> <span class="p">(</span><span class="n">ratio</span> <span class="o">&gt;=</span> <span class="mf">0.0001</span><span class="p">):</span>
                    <span class="c">## Successful iteration.  Update x, fvec, and their norms</span>
                    <span class="n">x</span> <span class="o">=</span> <span class="n">wa2</span>
                    <span class="n">wa2</span> <span class="o">=</span> <span class="n">diag</span> <span class="o">*</span> <span class="n">x</span>
                    <span class="n">fvec</span> <span class="o">=</span> <span class="n">wa4</span>
                    <span class="n">xnorm</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">enorm</span><span class="p">(</span><span class="n">wa2</span><span class="p">)</span>
                    <span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span> <span class="o">=</span> <span class="n">fnorm1</span>
                    <span class="bp">self</span><span class="o">.</span><span class="n">niter</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">niter</span> <span class="o">+</span> <span class="mi">1</span>
                
                <span class="c">## Tests for convergence</span>
                <span class="k">if</span> <span class="p">((</span><span class="nb">abs</span><span class="p">(</span><span class="n">actred</span><span class="p">)</span> <span class="o">&lt;=</span> <span class="n">ftol</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">prered</span> <span class="o">&lt;=</span> <span class="n">ftol</span><span class="p">)</span>
                     <span class="ow">and</span> <span class="p">(</span><span class="mf">0.5</span> <span class="o">*</span> <span class="n">ratio</span> <span class="o">&lt;=</span> <span class="mi">1</span><span class="p">)):</span> <span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">=</span> <span class="mi">1</span>
                <span class="k">if</span> <span class="n">delta</span> <span class="o">&lt;=</span> <span class="n">xtol</span><span class="o">*</span><span class="n">xnorm</span><span class="p">:</span> <span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">=</span> <span class="mi">2</span>
                <span class="k">if</span> <span class="p">((</span><span class="nb">abs</span><span class="p">(</span><span class="n">actred</span><span class="p">)</span> <span class="o">&lt;=</span> <span class="n">ftol</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">prered</span> <span class="o">&lt;=</span> <span class="n">ftol</span><span class="p">)</span>
                     <span class="ow">and</span> <span class="p">(</span><span class="mf">0.5</span> <span class="o">*</span> <span class="n">ratio</span> <span class="o">&lt;=</span> <span class="mi">1</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">==</span> <span class="mi">2</span><span class="p">)):</span> <span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">=</span> <span class="mi">3</span>
                <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">):</span> <span class="k">break</span>
                
                <span class="c">## Tests for termination and stringent tolerances</span>
                <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">niter</span> <span class="o">&gt;=</span> <span class="n">maxiter</span><span class="p">):</span> <span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">=</span> <span class="mi">5</span>
                <span class="k">if</span> <span class="p">((</span><span class="nb">abs</span><span class="p">(</span><span class="n">actred</span><span class="p">)</span> <span class="o">&lt;=</span> <span class="n">machep</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">prered</span> <span class="o">&lt;=</span> <span class="n">machep</span><span class="p">)</span>
                    <span class="ow">and</span> <span class="p">(</span><span class="mf">0.5</span><span class="o">*</span><span class="n">ratio</span> <span class="o">&lt;=</span> <span class="mi">1</span><span class="p">)):</span> <span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">=</span> <span class="mi">6</span>
                <span class="k">if</span> <span class="n">delta</span> <span class="o">&lt;=</span> <span class="n">machep</span><span class="o">*</span><span class="n">xnorm</span><span class="p">:</span> <span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">=</span> <span class="mi">7</span>
                <span class="k">if</span> <span class="n">gnorm</span> <span class="o">&lt;=</span> <span class="n">machep</span><span class="p">:</span> <span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">=</span> <span class="mi">8</span>
                <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">):</span> <span class="k">break</span>
                
                <span class="c">## End of inner loop. Repeat if iteration unsuccessful</span>
                <span class="k">if</span> <span class="p">(</span><span class="n">ratio</span> <span class="o">&gt;=</span> <span class="mf">0.0001</span><span class="p">):</span> <span class="k">break</span>

                <span class="c">## Check for over/underflow</span>
                <span class="k">if</span> <span class="o">~</span><span class="n">numpy</span><span class="o">.</span><span class="n">all</span><span class="p">(</span><span class="n">numpy</span><span class="o">.</span><span class="n">isfinite</span><span class="p">(</span><span class="n">wa1</span><span class="p">)</span> <span class="o">&amp;</span> <span class="n">numpy</span><span class="o">.</span><span class="n">isfinite</span><span class="p">(</span><span class="n">wa2</span><span class="p">)</span> <span class="o">&amp;</span>
                            <span class="n">numpy</span><span class="o">.</span><span class="n">isfinite</span><span class="p">(</span><span class="n">x</span><span class="p">))</span> <span class="ow">or</span> <span class="o">~</span><span class="n">numpy</span><span class="o">.</span><span class="n">isfinite</span><span class="p">(</span><span class="n">ratio</span><span class="p">):</span>
                    <span class="n">errmsg</span> <span class="o">=</span> <span class="p">(</span><span class="s">&#39;&#39;&#39;ERROR: parameter or function value(s) have become </span>
<span class="s">                        &#39;infinite; check model function for over- &#39;and underflow&#39;&#39;&#39;</span><span class="p">)</span>
                    <span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">=</span> <span class="o">-</span><span class="mi">16</span>
                    <span class="k">break</span>
                <span class="c">##wh = where(finite(wa1) EQ 0 OR finite(wa2) EQ 0 OR finite(x) EQ 0, ct)</span>
                <span class="c">##if ct GT 0 OR finite(ratio) EQ 0 then begin</span>

            <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">):</span> <span class="k">break</span><span class="p">;</span>
        <span class="c">## End of outer loop.</span>

        <span class="n">catch_msg</span> <span class="o">=</span> <span class="s">&#39;in the termination phase&#39;</span>
        <span class="c">## Termination, either normal or user imposed.</span>
        <span class="k">if</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="p">)</span> <span class="o">==</span> <span class="mi">0</span><span class="p">):</span>
            <span class="k">return</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">nfree</span> <span class="o">==</span> <span class="mi">0</span><span class="p">):</span> <span class="bp">self</span><span class="o">.</span><span class="n">params</span> <span class="o">=</span> <span class="n">xall</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
        <span class="k">else</span><span class="p">:</span> <span class="n">put</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="p">,</span> <span class="n">ifree</span><span class="p">,</span> <span class="n">x</span><span class="p">)</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">nprint</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span>
            <span class="n">catch_msg</span> <span class="o">=</span> <span class="s">&#39;calling &#39;</span> <span class="o">+</span> <span class="nb">str</span><span class="p">(</span><span class="n">fcn</span><span class="p">)</span>
            <span class="p">[</span><span class="n">status</span><span class="p">,</span> <span class="n">fvec</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">call</span><span class="p">(</span><span class="n">fcn</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="p">,</span> <span class="n">functkw</span><span class="p">)</span>
            <span class="n">catch_msg</span> <span class="o">=</span> <span class="s">&#39;in the termination phase&#39;</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">enorm</span><span class="p">(</span><span class="n">fvec</span><span class="p">)</span>

        <span class="k">if</span> <span class="p">((</span><span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">fnorm1</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">)):</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">max</span><span class="p">([</span><span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span><span class="p">,</span> <span class="n">fnorm1</span><span class="p">])</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">fnorm</span><span class="o">**</span><span class="mf">2.</span>

        <span class="bp">self</span><span class="o">.</span><span class="n">covar</span> <span class="o">=</span> <span class="bp">None</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">perror</span> <span class="o">=</span> <span class="bp">None</span>
        <span class="c">## (very carefully) set the covariance matrix COVAR</span>
        <span class="k">if</span> <span class="p">((</span><span class="bp">self</span><span class="o">.</span><span class="n">status</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">nocovar</span><span class="o">==</span><span class="mi">0</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">n</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">)</span>
                       <span class="ow">and</span> <span class="p">(</span><span class="n">fjac</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">ipvt</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">)):</span>
            <span class="n">sz</span> <span class="o">=</span> <span class="n">shape</span><span class="p">(</span><span class="n">fjac</span><span class="p">)</span>
            <span class="k">if</span> <span class="p">((</span><span class="n">n</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">sz</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">&gt;=</span> <span class="n">n</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">sz</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">&gt;=</span> <span class="n">n</span><span class="p">)</span>
                <span class="ow">and</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">ipvt</span><span class="p">)</span> <span class="o">&gt;=</span> <span class="n">n</span><span class="p">)):</span>

                <span class="n">catch_msg</span> <span class="o">=</span> <span class="s">&#39;computing the covariance matrix&#39;</span>
                <span class="n">cv</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">calc_covar</span><span class="p">(</span><span class="n">fjac</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">n</span><span class="p">,</span><span class="mi">0</span><span class="p">:</span><span class="n">n</span><span class="p">],</span> <span class="n">ipvt</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">n</span><span class="p">])</span>
                <span class="n">cv</span><span class="o">.</span><span class="n">shape</span> <span class="o">=</span> <span class="p">[</span><span class="n">n</span><span class="p">,</span> <span class="n">n</span><span class="p">]</span>
                <span class="n">nn</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">xall</span><span class="p">)</span>

                <span class="c">## Fill in actual covariance matrix, accounting for fixed</span>
                <span class="c">## parameters.</span>
                <span class="bp">self</span><span class="o">.</span><span class="n">covar</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">([</span><span class="n">nn</span><span class="p">,</span> <span class="n">nn</span><span class="p">],</span> <span class="n">dtype</span><span class="o">=</span><span class="nb">float</span><span class="p">)</span>
                <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
                    <span class="bp">self</span><span class="o">.</span><span class="n">covar</span><span class="p">[</span><span class="n">ifree</span><span class="p">,</span><span class="n">ifree</span><span class="p">[</span><span class="n">i</span><span class="p">]]</span> <span class="o">=</span> <span class="n">cv</span><span class="p">[:,</span><span class="n">i</span><span class="p">]</span>

                <span class="c">## Compute errors in parameters</span>
                <span class="n">catch_msg</span> <span class="o">=</span> <span class="s">&#39;computing parameter errors&#39;</span>
                <span class="bp">self</span><span class="o">.</span><span class="n">perror</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">nn</span><span class="p">,</span> <span class="n">dtype</span><span class="o">=</span><span class="nb">float</span><span class="p">)</span>
                <span class="n">d</span> <span class="o">=</span> <span class="n">diagonal</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">covar</span><span class="p">)</span>
                <span class="n">wh</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">(</span><span class="n">d</span> <span class="o">&gt;=</span> <span class="mi">0</span><span class="p">))[</span><span class="mi">0</span><span class="p">]</span>
                <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">wh</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
                    <span class="n">put</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">perror</span><span class="p">,</span> <span class="n">wh</span><span class="p">,</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">take</span><span class="p">(</span><span class="n">d</span><span class="p">,</span> <span class="n">wh</span><span class="p">)))</span>
        <span class="k">return</span>


    <span class="k">def</span> <span class="nf">__str__</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
        <span class="k">return</span> <span class="p">{</span><span class="s">&#39;params&#39;</span><span class="p">:</span> <span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="p">,</span>
               <span class="s">&#39;niter&#39;</span><span class="p">:</span> <span class="bp">self</span><span class="o">.</span><span class="n">niter</span><span class="p">,</span>
               <span class="s">&#39;params&#39;</span><span class="p">:</span> <span class="bp">self</span><span class="o">.</span><span class="n">params</span><span class="p">,</span>
               <span class="s">&#39;covar&#39;</span><span class="p">:</span> <span class="bp">self</span><span class="o">.</span><span class="n">covar</span><span class="p">,</span>
               <span class="s">&#39;perror&#39;</span><span class="p">:</span> <span class="bp">self</span><span class="o">.</span><span class="n">perror</span><span class="p">,</span>
               <span class="s">&#39;status&#39;</span><span class="p">:</span> <span class="bp">self</span><span class="o">.</span><span class="n">status</span><span class="p">,</span>
               <span class="s">&#39;debug&#39;</span><span class="p">:</span> <span class="bp">self</span><span class="o">.</span><span class="n">debug</span><span class="p">,</span>
               <span class="s">&#39;errmsg&#39;</span><span class="p">:</span> <span class="bp">self</span><span class="o">.</span><span class="n">errmsg</span><span class="p">,</span>
               <span class="s">&#39;fastnorm&#39;</span><span class="p">:</span> <span class="bp">self</span><span class="o">.</span><span class="n">fastnorm</span><span class="p">,</span>
               <span class="s">&#39;nfev&#39;</span><span class="p">:</span> <span class="bp">self</span><span class="o">.</span><span class="n">nfev</span><span class="p">,</span>
               <span class="s">&#39;damp&#39;</span><span class="p">:</span> <span class="bp">self</span><span class="o">.</span><span class="n">damp</span>
               <span class="c">#,&#39;machar&#39;:self.machar</span>
               <span class="p">}</span><span class="o">.</span><span class="n">__str__</span><span class="p">()</span>

    <span class="c">## Default procedure to be called every iteration.  It simply prints</span>
    <span class="c">## the parameter values.</span>
<div class="viewcode-block" id="mpfit.defiter"><a class="viewcode-back" href="../../agpy.html#agpy.mpfit_custom.mpfit.defiter">[docs]</a>    <span class="k">def</span> <span class="nf">defiter</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">fcn</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="nb">iter</span><span class="p">,</span> <span class="n">fnorm</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">functkw</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span>
                       <span class="n">quiet</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">iterstop</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">parinfo</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span>
                       <span class="n">format</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">pformat</span><span class="o">=</span><span class="s">&#39;</span><span class="si">%.10g</span><span class="s">&#39;</span><span class="p">,</span> <span class="n">dof</span><span class="o">=</span><span class="mi">1</span><span class="p">):</span>

        <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">debug</span><span class="p">):</span> <span class="k">print</span> <span class="s">&#39;Entering defiter...&#39;</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">quiet</span><span class="p">):</span> <span class="k">return</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">fnorm</span> <span class="o">==</span> <span class="bp">None</span><span class="p">):</span>
            <span class="p">[</span><span class="n">status</span><span class="p">,</span> <span class="n">fvec</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">call</span><span class="p">(</span><span class="n">fcn</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">functkw</span><span class="p">)</span>
            <span class="n">fnorm</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">enorm</span><span class="p">(</span><span class="n">fvec</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span>

        <span class="c">## Determine which parameters to print</span>
        <span class="n">nprint</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
        <span class="k">print</span> <span class="s">&quot;Iter &quot;</span><span class="p">,</span> <span class="p">(</span><span class="s">&#39;</span><span class="si">%6i</span><span class="s">&#39;</span> <span class="o">%</span> <span class="nb">iter</span><span class="p">),</span><span class="s">&quot;   CHI-SQUARE = &quot;</span><span class="p">,(</span><span class="s">&#39;</span><span class="si">%.10g</span><span class="s">&#39;</span> <span class="o">%</span> <span class="n">fnorm</span><span class="p">),</span><span class="s">&quot; DOF = &quot;</span><span class="p">,</span> <span class="p">(</span><span class="s">&#39;</span><span class="si">%i</span><span class="s">&#39;</span> <span class="o">%</span> <span class="n">dof</span><span class="p">)</span>
        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nprint</span><span class="p">):</span>
            <span class="k">if</span> <span class="p">(</span><span class="n">parinfo</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">parinfo</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">has_key</span><span class="p">(</span><span class="s">&#39;parname&#39;</span><span class="p">)):</span>
                <span class="n">p</span> <span class="o">=</span> <span class="s">&#39;   &#39;</span> <span class="o">+</span> <span class="n">parinfo</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="s">&#39;parname&#39;</span><span class="p">]</span> <span class="o">+</span> <span class="s">&#39; = &#39;</span>
            <span class="k">else</span><span class="p">:</span>
                <span class="n">p</span> <span class="o">=</span> <span class="s">&#39;   P&#39;</span> <span class="o">+</span> <span class="nb">str</span><span class="p">(</span><span class="n">i</span><span class="p">)</span> <span class="o">+</span> <span class="s">&#39; = &#39;</span>
            <span class="k">if</span> <span class="p">(</span><span class="n">parinfo</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">parinfo</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">has_key</span><span class="p">(</span><span class="s">&#39;mpprint&#39;</span><span class="p">)):</span>
                <span class="n">iprint</span> <span class="o">=</span> <span class="n">parinfo</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="s">&#39;mpprint&#39;</span><span class="p">]</span>
            <span class="k">else</span><span class="p">:</span>
                <span class="n">iprint</span> <span class="o">=</span> <span class="mi">1</span>
            <span class="k">if</span> <span class="p">(</span><span class="n">iprint</span><span class="p">):</span>
                <span class="k">print</span> <span class="n">p</span> <span class="o">+</span> <span class="p">(</span><span class="n">pformat</span> <span class="o">%</span> <span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">])</span> <span class="o">+</span> <span class="s">&#39;  &#39;</span>
        <span class="k">return</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>

    <span class="c">##  DO_ITERSTOP:</span>
    <span class="c">##  if keyword_set(iterstop) then begin</span>
    <span class="c">##    k = get_kbrd(0)</span>
    <span class="c">##    if k EQ string(byte(7)) then begin</span>
    <span class="c">##        message, &#39;WARNING: minimization not complete&#39;, /info</span>
    <span class="c">##        print, &#39;Do you want to terminate this procedure? (y/n)&#39;, $</span>
    <span class="c">##          format=&#39;(A,$)&#39;</span>
    <span class="c">##        k = &#39;&#39;</span>
    <span class="c">##        read, k</span>
    <span class="c">##        if strupcase(strmid(k,0,1)) EQ &#39;Y&#39; then begin</span>
    <span class="c">##            message, &#39;WARNING: Procedure is terminating.&#39;, /info</span>
    <span class="c">##            mperr = -1</span>
    <span class="c">##        endif</span>
    <span class="c">##    endif</span>
    <span class="c">##  endif</span>
    
    
    <span class="c">## Procedure to parse the parameter values in PARINFO, which is a list of dictionaries</span></div>
<div class="viewcode-block" id="mpfit.parinfo"><a class="viewcode-back" href="../../agpy.html#agpy.mpfit_custom.mpfit.parinfo">[docs]</a>    <span class="k">def</span> <span class="nf">parinfo</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">parinfo</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">key</span><span class="o">=</span><span class="s">&#39;a&#39;</span><span class="p">,</span> <span class="n">default</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">n</span><span class="o">=</span><span class="mi">0</span><span class="p">):</span>
        <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">debug</span><span class="p">):</span> <span class="k">print</span> <span class="s">&#39;Entering parinfo...&#39;</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">n</span> <span class="o">==</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">parinfo</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">):</span> <span class="n">n</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">parinfo</span><span class="p">)</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">n</span> <span class="o">==</span> <span class="mi">0</span><span class="p">):</span>
            <span class="n">values</span> <span class="o">=</span> <span class="n">default</span>
            <span class="k">return</span><span class="p">(</span><span class="n">values</span><span class="p">)</span>
        <span class="n">values</span> <span class="o">=</span> <span class="p">[]</span>
        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
            <span class="k">if</span> <span class="p">((</span><span class="n">parinfo</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">parinfo</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">has_key</span><span class="p">(</span><span class="n">key</span><span class="p">))):</span>
                <span class="n">values</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">parinfo</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">key</span><span class="p">])</span>
            <span class="k">else</span><span class="p">:</span>
                <span class="n">values</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">default</span><span class="p">)</span>

        <span class="c"># Convert to numeric arrays if possible</span>
        <span class="n">test</span> <span class="o">=</span> <span class="n">default</span>
        <span class="k">if</span> <span class="p">(</span><span class="nb">type</span><span class="p">(</span><span class="n">default</span><span class="p">)</span> <span class="o">==</span> <span class="n">types</span><span class="o">.</span><span class="n">ListType</span><span class="p">):</span> <span class="n">test</span><span class="o">=</span><span class="n">default</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
        <span class="k">if</span> <span class="p">(</span><span class="nb">isinstance</span><span class="p">(</span><span class="n">test</span><span class="p">,</span> <span class="n">types</span><span class="o">.</span><span class="n">IntType</span><span class="p">)):</span>
            <span class="n">values</span> <span class="o">=</span> <span class="n">asarray</span><span class="p">(</span><span class="n">values</span><span class="p">,</span> <span class="nb">int</span><span class="p">)</span>
        <span class="k">elif</span> <span class="p">(</span><span class="nb">isinstance</span><span class="p">(</span><span class="n">test</span><span class="p">,</span> <span class="n">types</span><span class="o">.</span><span class="n">FloatType</span><span class="p">)):</span>
            <span class="n">values</span> <span class="o">=</span> <span class="n">asarray</span><span class="p">(</span><span class="n">values</span><span class="p">,</span> <span class="nb">float</span><span class="p">)</span>
        <span class="k">return</span><span class="p">(</span><span class="n">values</span><span class="p">)</span>
    
    <span class="c">## Call user function or procedure, with _EXTRA or not, with</span>
    <span class="c">## derivatives or not.</span></div>
<div class="viewcode-block" id="mpfit.call"><a class="viewcode-back" href="../../agpy.html#agpy.mpfit_custom.mpfit.call">[docs]</a>    <span class="k">def</span> <span class="nf">call</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">fcn</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">functkw</span><span class="p">,</span> <span class="n">fjac</span><span class="o">=</span><span class="bp">None</span><span class="p">):</span>
        <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">debug</span><span class="p">):</span> <span class="k">print</span> <span class="s">&#39;Entering call...&#39;</span>
        <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">qanytied</span><span class="p">):</span> <span class="n">x</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">tie</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">ptied</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">nfev</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nfev</span> <span class="o">+</span> <span class="mi">1</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">fjac</span> <span class="o">==</span> <span class="bp">None</span><span class="p">):</span>
            <span class="p">[</span><span class="n">status</span><span class="p">,</span> <span class="n">f</span><span class="p">]</span> <span class="o">=</span> <span class="n">fcn</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">fjac</span><span class="o">=</span><span class="n">fjac</span><span class="p">,</span> <span class="o">**</span><span class="n">functkw</span><span class="p">)</span>
            <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">damp</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span>
                <span class="c">## Apply the damping if requested.  This replaces the residuals</span>
                <span class="c">## with their hyperbolic tangent.  Thus residuals larger than</span>
                <span class="c">## DAMP are essentially clipped.</span>
                <span class="n">f</span> <span class="o">=</span> <span class="n">tanh</span><span class="p">(</span><span class="n">f</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">damp</span><span class="p">)</span>
            <span class="k">return</span><span class="p">([</span><span class="n">status</span><span class="p">,</span> <span class="n">f</span><span class="p">])</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="k">return</span><span class="p">(</span><span class="n">fcn</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">fjac</span><span class="o">=</span><span class="n">fjac</span><span class="p">,</span> <span class="o">**</span><span class="n">functkw</span><span class="p">))</span>
    
    </div>
<div class="viewcode-block" id="mpfit.enorm"><a class="viewcode-back" href="../../agpy.html#agpy.mpfit_custom.mpfit.enorm">[docs]</a>    <span class="k">def</span> <span class="nf">enorm</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">vec</span><span class="p">):</span>

        <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">debug</span><span class="p">):</span> <span class="k">print</span> <span class="s">&#39;Entering enorm...&#39;</span>
        <span class="c">## NOTE: it turns out that, for systems that have a lot of data</span>
        <span class="c">## points, this routine is a big computing bottleneck.  The extended</span>
        <span class="c">## computations that need to be done cannot be effectively</span>
        <span class="c">## vectorized.  The introduction of the FASTNORM configuration</span>
        <span class="c">## parameter allows the user to select a faster routine, which is</span>
        <span class="c">## based on TOTAL() alone.</span>

        <span class="c"># Very simple-minded sum-of-squares</span>
        <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">fastnorm</span><span class="p">):</span>
            <span class="n">ans</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="nb">sum</span><span class="p">(</span><span class="n">vec</span><span class="o">*</span><span class="n">vec</span><span class="p">))</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="n">agiant</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">machar</span><span class="o">.</span><span class="n">rgiant</span> <span class="o">/</span> <span class="nb">len</span><span class="p">(</span><span class="n">vec</span><span class="p">)</span>
            <span class="n">adwarf</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">machar</span><span class="o">.</span><span class="n">rdwarf</span> <span class="o">*</span> <span class="nb">len</span><span class="p">(</span><span class="n">vec</span><span class="p">)</span>

            <span class="c">## This is hopefully a compromise between speed and robustness.</span>
            <span class="c">## Need to do this because of the possibility of over- or underflow.</span>
            <span class="n">mx</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">max</span><span class="p">(</span><span class="n">vec</span><span class="p">)</span>
            <span class="n">mn</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">min</span><span class="p">(</span><span class="n">vec</span><span class="p">)</span>
            <span class="n">mx</span> <span class="o">=</span> <span class="nb">max</span><span class="p">([</span><span class="nb">abs</span><span class="p">(</span><span class="n">mx</span><span class="p">),</span> <span class="nb">abs</span><span class="p">(</span><span class="n">mn</span><span class="p">)])</span>
            <span class="k">if</span> <span class="n">mx</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span> <span class="k">return</span><span class="p">(</span><span class="n">vec</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">*</span><span class="mf">0.</span><span class="p">)</span>
            <span class="k">if</span> <span class="n">mx</span> <span class="o">&gt;</span> <span class="n">agiant</span> <span class="ow">or</span> <span class="n">mx</span> <span class="o">&lt;</span> <span class="n">adwarf</span><span class="p">:</span>
                <span class="n">ans</span> <span class="o">=</span> <span class="n">mx</span> <span class="o">*</span> <span class="n">sqrt</span><span class="p">(</span><span class="nb">sum</span><span class="p">((</span><span class="n">vec</span><span class="o">/</span><span class="n">mx</span><span class="p">)</span><span class="o">*</span><span class="p">(</span><span class="n">vec</span><span class="o">/</span><span class="n">mx</span><span class="p">)))</span>
            <span class="k">else</span><span class="p">:</span>
                <span class="n">ans</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="nb">sum</span><span class="p">(</span><span class="n">vec</span><span class="o">*</span><span class="n">vec</span><span class="p">))</span>

        <span class="k">return</span><span class="p">(</span><span class="n">ans</span><span class="p">)</span>
    
    </div>
<div class="viewcode-block" id="mpfit.fdjac2"><a class="viewcode-back" href="../../agpy.html#agpy.mpfit_custom.mpfit.fdjac2">[docs]</a>    <span class="k">def</span> <span class="nf">fdjac2</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">fcn</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">fvec</span><span class="p">,</span> <span class="n">step</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">ulimited</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">ulimit</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">dside</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span>
               <span class="n">epsfcn</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">autoderivative</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span>
               <span class="n">functkw</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">xall</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">ifree</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">dstep</span><span class="o">=</span><span class="bp">None</span><span class="p">):</span>

        <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">debug</span><span class="p">):</span> <span class="k">print</span> <span class="s">&#39;Entering fdjac2...&#39;</span>
        <span class="n">machep</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">machar</span><span class="o">.</span><span class="n">machep</span>
        <span class="k">if</span> <span class="n">epsfcn</span> <span class="o">==</span> <span class="bp">None</span><span class="p">:</span>  <span class="n">epsfcn</span> <span class="o">=</span> <span class="n">machep</span>
        <span class="k">if</span> <span class="n">xall</span> <span class="o">==</span> <span class="bp">None</span><span class="p">:</span>    <span class="n">xall</span> <span class="o">=</span> <span class="n">x</span>
        <span class="k">if</span> <span class="n">ifree</span> <span class="o">==</span> <span class="bp">None</span><span class="p">:</span>   <span class="n">ifree</span> <span class="o">=</span> <span class="n">arange</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">xall</span><span class="p">))</span>
        <span class="k">if</span> <span class="n">step</span> <span class="o">==</span> <span class="bp">None</span><span class="p">:</span>    <span class="n">step</span> <span class="o">=</span> <span class="n">x</span> <span class="o">*</span> <span class="mf">0.</span>
        <span class="n">nall</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">xall</span><span class="p">)</span>

        <span class="n">eps</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">numpy</span><span class="o">.</span><span class="n">max</span><span class="p">([</span><span class="n">epsfcn</span><span class="p">,</span> <span class="n">machep</span><span class="p">]))</span>
        <span class="n">m</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">fvec</span><span class="p">)</span>
        <span class="n">n</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>

        <span class="c">## Compute analytical derivative if requested</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">autoderivative</span> <span class="o">==</span> <span class="mi">0</span><span class="p">):</span>
            <span class="n">mperr</span> <span class="o">=</span> <span class="mi">0</span>
            <span class="n">fjac</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">nall</span><span class="p">,</span> <span class="n">dtype</span><span class="o">=</span><span class="nb">float</span><span class="p">)</span>
            <span class="n">Put</span><span class="p">(</span><span class="n">fjac</span><span class="p">,</span> <span class="n">ifree</span><span class="p">,</span> <span class="mf">1.0</span><span class="p">)</span>  <span class="c">## Specify which parameters need derivatives</span>
            <span class="p">[</span><span class="n">status</span><span class="p">,</span> <span class="n">fp</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">call</span><span class="p">(</span><span class="n">fcn</span><span class="p">,</span> <span class="n">xall</span><span class="p">,</span> <span class="n">functkw</span><span class="p">,</span> <span class="n">fjac</span><span class="o">=</span><span class="n">fjac</span><span class="p">)</span>

            <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">fjac</span><span class="p">)</span> <span class="o">!=</span> <span class="n">m</span><span class="o">*</span><span class="n">nall</span><span class="p">:</span>
                <span class="k">print</span> <span class="s">&#39;ERROR: Derivative matrix was not computed properly.&#39;</span>
                <span class="k">return</span><span class="p">(</span><span class="bp">None</span><span class="p">)</span>

            <span class="c">## This definition is consistent with CURVEFIT</span>
            <span class="c">## Sign error found (thanks Jesus Fernandez &lt;fernande@irm.chu-caen.fr&gt;)</span>
            <span class="n">fjac</span><span class="o">.</span><span class="n">shape</span> <span class="o">=</span> <span class="p">[</span><span class="n">m</span><span class="p">,</span><span class="n">nall</span><span class="p">]</span>
            <span class="n">fjac</span> <span class="o">=</span> <span class="o">-</span><span class="n">fjac</span>

            <span class="c">## Select only the free parameters</span>
            <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">ifree</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">nall</span><span class="p">:</span>
                <span class="n">fjac</span> <span class="o">=</span> <span class="n">fjac</span><span class="p">[:,</span><span class="n">ifree</span><span class="p">]</span>
                <span class="n">fjac</span><span class="o">.</span><span class="n">shape</span> <span class="o">=</span> <span class="p">[</span><span class="n">m</span><span class="p">,</span> <span class="n">n</span><span class="p">]</span>
                <span class="k">return</span><span class="p">(</span><span class="n">fjac</span><span class="p">)</span>

        <span class="n">fjac</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">([</span><span class="n">m</span><span class="p">,</span> <span class="n">n</span><span class="p">],</span> <span class="n">dtype</span><span class="o">=</span><span class="nb">float</span><span class="p">)</span>

        <span class="n">h</span> <span class="o">=</span> <span class="n">eps</span> <span class="o">*</span> <span class="nb">abs</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>

        <span class="c">## if STEP is given, use that</span>
        <span class="k">if</span> <span class="n">step</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">:</span>
            <span class="n">stepi</span> <span class="o">=</span> <span class="n">take</span><span class="p">(</span><span class="n">step</span><span class="p">,</span> <span class="n">ifree</span><span class="p">)</span>
            <span class="n">wh</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">(</span><span class="n">stepi</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">))[</span><span class="mi">0</span><span class="p">]</span>
            <span class="k">if</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">wh</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span> <span class="n">put</span><span class="p">(</span><span class="n">h</span><span class="p">,</span> <span class="n">wh</span><span class="p">,</span> <span class="n">take</span><span class="p">(</span><span class="n">stepi</span><span class="p">,</span> <span class="n">wh</span><span class="p">))</span>

        <span class="c">## if relative step is given, use that</span>
        <span class="k">if</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">dstep</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span>
            <span class="n">dstepi</span> <span class="o">=</span> <span class="n">take</span><span class="p">(</span><span class="n">dstep</span><span class="p">,</span> <span class="n">ifree</span><span class="p">)</span>
            <span class="n">wh</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">(</span><span class="n">dstepi</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">))[</span><span class="mi">0</span><span class="p">]</span>
            <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">wh</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span> <span class="n">put</span><span class="p">(</span><span class="n">h</span><span class="p">,</span> <span class="n">wh</span><span class="p">,</span> <span class="nb">abs</span><span class="p">(</span><span class="n">take</span><span class="p">(</span><span class="n">dstepi</span><span class="p">,</span><span class="n">wh</span><span class="p">)</span><span class="o">*</span><span class="n">take</span><span class="p">(</span><span class="n">x</span><span class="p">,</span><span class="n">wh</span><span class="p">)))</span>

        <span class="c">## In case any of the step values are zero</span>
        <span class="n">wh</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">(</span><span class="n">h</span> <span class="o">==</span> <span class="mi">0</span><span class="p">))[</span><span class="mi">0</span><span class="p">]</span>
        <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">wh</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span> <span class="n">put</span><span class="p">(</span><span class="n">h</span><span class="p">,</span> <span class="n">wh</span><span class="p">,</span> <span class="n">eps</span><span class="p">)</span>

        <span class="c">## Reverse the sign of the step if we are up against the parameter</span>
        <span class="c">## limit, or if the user requested it.</span>
        <span class="n">mask</span> <span class="o">=</span> <span class="n">dside</span><span class="p">[</span><span class="n">ifree</span><span class="p">]</span> <span class="o">==</span> <span class="o">-</span><span class="mi">1</span>
        <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">ulimited</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span> <span class="ow">and</span> <span class="nb">len</span><span class="p">(</span><span class="n">ulimit</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
            <span class="n">mask</span> <span class="o">=</span> <span class="n">logical_or</span><span class="p">((</span><span class="n">mask</span><span class="o">!=</span><span class="mi">0</span><span class="p">),</span> <span class="n">logical_and</span><span class="p">((</span><span class="n">ulimited</span><span class="o">!=</span><span class="mi">0</span><span class="p">),(</span><span class="n">x</span> <span class="o">&gt;</span> <span class="n">ulimit</span><span class="o">-</span><span class="n">h</span><span class="p">)))</span>
            <span class="n">wh</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">(</span><span class="n">mask</span><span class="p">))[</span><span class="mi">0</span><span class="p">]</span>
            <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">wh</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span> <span class="n">put</span><span class="p">(</span><span class="n">h</span><span class="p">,</span> <span class="n">wh</span><span class="p">,</span> <span class="o">-</span><span class="n">take</span><span class="p">(</span><span class="n">h</span><span class="p">,</span> <span class="n">wh</span><span class="p">))</span>
        <span class="c">## Loop through parameters, computing the derivative for each</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
            <span class="n">xp</span> <span class="o">=</span> <span class="n">xall</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
            <span class="n">xp</span><span class="p">[</span><span class="n">ifree</span><span class="p">[</span><span class="n">j</span><span class="p">]]</span> <span class="o">=</span> <span class="n">xp</span><span class="p">[</span><span class="n">ifree</span><span class="p">[</span><span class="n">j</span><span class="p">]]</span> <span class="o">+</span> <span class="n">h</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
            <span class="p">[</span><span class="n">status</span><span class="p">,</span> <span class="n">fp</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">call</span><span class="p">(</span><span class="n">fcn</span><span class="p">,</span> <span class="n">xp</span><span class="p">,</span> <span class="n">functkw</span><span class="p">)</span>
            <span class="k">if</span> <span class="p">(</span><span class="n">status</span> <span class="o">&lt;</span> <span class="mi">0</span><span class="p">):</span> <span class="k">return</span><span class="p">(</span><span class="bp">None</span><span class="p">)</span>

            <span class="k">if</span> <span class="nb">abs</span><span class="p">(</span><span class="n">dside</span><span class="p">[</span><span class="n">j</span><span class="p">])</span> <span class="o">&lt;=</span> <span class="mi">1</span><span class="p">:</span>
                <span class="c">## COMPUTE THE ONE-SIDED DERIVATIVE</span>
                <span class="c">## Note optimization fjac(0:*,j)</span>
                <span class="n">fjac</span><span class="p">[</span><span class="mi">0</span><span class="p">:,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">fp</span><span class="o">-</span><span class="n">fvec</span><span class="p">)</span><span class="o">/</span><span class="n">h</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>

            <span class="k">else</span><span class="p">:</span>
                <span class="c">## COMPUTE THE TWO-SIDED DERIVATIVE</span>
                <span class="n">xp</span><span class="p">[</span><span class="n">ifree</span><span class="p">[</span><span class="n">j</span><span class="p">]]</span> <span class="o">=</span> <span class="n">xall</span><span class="p">[</span><span class="n">ifree</span><span class="p">[</span><span class="n">j</span><span class="p">]]</span> <span class="o">-</span> <span class="n">h</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>

                <span class="n">mperr</span> <span class="o">=</span> <span class="mi">0</span>
                <span class="p">[</span><span class="n">status</span><span class="p">,</span> <span class="n">fm</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">call</span><span class="p">(</span><span class="n">fcn</span><span class="p">,</span> <span class="n">xp</span><span class="p">,</span> <span class="n">functkw</span><span class="p">)</span>
                <span class="k">if</span> <span class="p">(</span><span class="n">status</span> <span class="o">&lt;</span> <span class="mi">0</span><span class="p">):</span> <span class="k">return</span><span class="p">(</span><span class="bp">None</span><span class="p">)</span>

                <span class="c">## Note optimization fjac(0:*,j)</span>
                <span class="n">fjac</span><span class="p">[</span><span class="mi">0</span><span class="p">:,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">fp</span><span class="o">-</span><span class="n">fm</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">h</span><span class="p">[</span><span class="n">j</span><span class="p">])</span>
        <span class="k">return</span><span class="p">(</span><span class="n">fjac</span><span class="p">)</span>
    
    
    
    <span class="c">#    Original FORTRAN documentation</span>
    <span class="c">#    **********</span>
    <span class="c">#</span>
    <span class="c">#    subroutine qrfac</span>
    <span class="c">#</span>
    <span class="c">#    this subroutine uses householder transformations with column</span>
    <span class="c">#    pivoting (optional) to compute a qr factorization of the</span>
    <span class="c">#    m by n matrix a. that is, qrfac determines an orthogonal</span>
    <span class="c">#    matrix q, a permutation matrix p, and an upper trapezoidal</span>
    <span class="c">#    matrix r with diagonal elements of nonincreasing magnitude,</span>
    <span class="c">#    such that a*p = q*r. the householder transformation for</span>
    <span class="c">#    column k, k = 1,2,...,min(m,n), is of the form</span>
    <span class="c">#</span>
    <span class="c">#                       t</span>
    <span class="c">#       i - (1/u(k))*u*u</span>
    <span class="c">#</span>
    <span class="c">#    where u has zeros in the first k-1 positions. the form of</span>
    <span class="c">#    this transformation and the method of pivoting first</span>
    <span class="c">#    appeared in the corresponding linpack subroutine.</span>
    <span class="c">#</span>
    <span class="c">#    the subroutine statement is</span>
    <span class="c">#</span>
    <span class="c">#   subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)</span>
    <span class="c">#</span>
    <span class="c">#    where</span>
    <span class="c">#</span>
    <span class="c">#   m is a positive integer input variable set to the number</span>
    <span class="c">#     of rows of a.</span>
    <span class="c">#</span>
    <span class="c">#   n is a positive integer input variable set to the number</span>
    <span class="c">#     of columns of a.</span>
    <span class="c">#</span>
    <span class="c">#   a is an m by n array. on input a contains the matrix for</span>
    <span class="c">#     which the qr factorization is to be computed. on output</span>
    <span class="c">#     the strict upper trapezoidal part of a contains the strict</span>
    <span class="c">#     upper trapezoidal part of r, and the lower trapezoidal</span>
    <span class="c">#     part of a contains a factored form of q (the non-trivial</span>
    <span class="c">#     elements of the u vectors described above).</span>
    <span class="c">#</span>
    <span class="c">#   lda is a positive integer input variable not less than m</span>
    <span class="c">#     which specifies the leading dimension of the array a.</span>
    <span class="c">#</span>
    <span class="c">#   pivot is a logical input variable. if pivot is set true,</span>
    <span class="c">#     then column pivoting is enforced. if pivot is set false,</span>
    <span class="c">#     then no column pivoting is done.</span>
    <span class="c">#</span>
    <span class="c">#   ipvt is an integer output array of length lipvt. ipvt</span>
    <span class="c">#     defines the permutation matrix p such that a*p = q*r.</span>
    <span class="c">#     column j of p is column ipvt(j) of the identity matrix.</span>
    <span class="c">#     if pivot is false, ipvt is not referenced.</span>
    <span class="c">#</span>
    <span class="c">#   lipvt is a positive integer input variable. if pivot is false,</span>
    <span class="c">#     then lipvt may be as small as 1. if pivot is true, then</span>
    <span class="c">#     lipvt must be at least n.</span>
    <span class="c">#</span>
    <span class="c">#   rdiag is an output array of length n which contains the</span>
    <span class="c">#     diagonal elements of r.</span>
    <span class="c">#</span>
    <span class="c">#   acnorm is an output array of length n which contains the</span>
    <span class="c">#     norms of the corresponding columns of the input matrix a.</span>
    <span class="c">#     if this information is not needed, then acnorm can coincide</span>
    <span class="c">#     with rdiag.</span>
    <span class="c">#</span>
    <span class="c">#   wa is a work array of length n. if pivot is false, then wa</span>
    <span class="c">#     can coincide with rdiag.</span>
    <span class="c">#</span>
    <span class="c">#    subprograms called</span>
    <span class="c">#</span>
    <span class="c">#   minpack-supplied ... dpmpar,enorm</span>
    <span class="c">#</span>
    <span class="c">#   fortran-supplied ... dmax1,dsqrt,min0</span>
    <span class="c">#</span>
    <span class="c">#    argonne national laboratory. minpack project. march 1980.</span>
    <span class="c">#    burton s. garbow, kenneth e. hillstrom, jorge j. more</span>
    <span class="c">#</span>
    <span class="c">#    **********</span>

    <span class="c"># NOTE: in IDL the factors appear slightly differently than described</span>
    <span class="c"># above.  The matrix A is still m x n where m &gt;= n.</span>
    <span class="c">#</span>
    <span class="c"># The &quot;upper&quot; triangular matrix R is actually stored in the strict</span>
    <span class="c"># lower left triangle of A under the standard notation of IDL.</span>
    <span class="c">#</span>
    <span class="c"># The reflectors that generate Q are in the upper trapezoid of A upon</span>
    <span class="c"># output.</span>
    <span class="c">#</span>
    <span class="c">#  EXAMPLE:  decompose the matrix [[9.,2.,6.],[4.,8.,7.]]</span>
    <span class="c">#   aa = [[9.,2.,6.],[4.,8.,7.]]</span>
    <span class="c">#   mpfit_qrfac, aa, aapvt, rdiag, aanorm</span>
    <span class="c">#    IDL&gt; print, aa</span>
    <span class="c">#         1.81818*   0.181818*   0.545455*</span>
    <span class="c">#        -8.54545+    1.90160*   0.432573*</span>
    <span class="c">#    IDL&gt; print, rdiag</span>
    <span class="c">#        -11.0000+   -7.48166+</span>
    <span class="c">#</span>
    <span class="c"># The components marked with a * are the components of the</span>
    <span class="c"># reflectors, and those marked with a + are components of R.</span>
    <span class="c">#</span>
    <span class="c"># To reconstruct Q and R we proceed as follows.  First R.</span>
    <span class="c">#   r = fltarr(m, n)</span>
    <span class="c">#   for i = 0, n-1 do r(0:i,i) = aa(0:i,i)  # fill in lower diag</span>
    <span class="c">#   r(lindgen(n)*(m+1)) = rdiag</span>
    <span class="c">#</span>
    <span class="c"># Next, Q, which are composed from the reflectors.  Each reflector v</span>
    <span class="c"># is taken from the upper trapezoid of aa, and converted to a matrix</span>
    <span class="c"># via (I - 2 vT . v / (v . vT)).</span>
    <span class="c">#</span>
    <span class="c">#   hh = ident                                  ## identity matrix</span>
    <span class="c">#   for i = 0, n-1 do begin</span>
    <span class="c">#   v = aa(*,i) &amp; if i GT 0 then v(0:i-1) = 0   ## extract reflector</span>
    <span class="c">#   hh = hh ## (ident - 2*(v # v)/total(v * v))  ## generate matrix</span>
    <span class="c">#   endfor</span>
    <span class="c">#</span>
    <span class="c"># Test the result:</span>
    <span class="c">#   IDL&gt; print, hh ## transpose(r)</span>
    <span class="c">#         9.00000     4.00000</span>
    <span class="c">#         2.00000     8.00000</span>
    <span class="c">#         6.00000     7.00000</span>
    <span class="c">#</span>
    <span class="c"># Note that it is usually never necessary to form the Q matrix</span>
    <span class="c"># explicitly, and MPFIT does not.</span>
    
</div>
<div class="viewcode-block" id="mpfit.qrfac"><a class="viewcode-back" href="../../agpy.html#agpy.mpfit_custom.mpfit.qrfac">[docs]</a>    <span class="k">def</span> <span class="nf">qrfac</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">a</span><span class="p">,</span> <span class="n">pivot</span><span class="o">=</span><span class="mi">0</span><span class="p">):</span>

        <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">debug</span><span class="p">):</span> <span class="k">print</span> <span class="s">&#39;Entering qrfac...&#39;</span>
        <span class="n">machep</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">machar</span><span class="o">.</span><span class="n">machep</span>
        <span class="n">sz</span> <span class="o">=</span> <span class="n">shape</span><span class="p">(</span><span class="n">a</span><span class="p">)</span>
        <span class="n">m</span> <span class="o">=</span> <span class="n">sz</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
        <span class="n">n</span> <span class="o">=</span> <span class="n">sz</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>

        <span class="c">## Compute the initial column norms and initialize arrays</span>
        <span class="n">acnorm</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">,</span> <span class="n">dtype</span><span class="o">=</span><span class="nb">float</span><span class="p">)</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
            <span class="n">acnorm</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">enorm</span><span class="p">(</span><span class="n">a</span><span class="p">[:,</span><span class="n">j</span><span class="p">])</span>
        <span class="n">rdiag</span> <span class="o">=</span> <span class="n">acnorm</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
        <span class="n">wa</span> <span class="o">=</span> <span class="n">rdiag</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
        <span class="n">ipvt</span> <span class="o">=</span> <span class="n">arange</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>

        <span class="c">## Reduce a to r with householder transformations</span>
        <span class="n">minmn</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">min</span><span class="p">([</span><span class="n">m</span><span class="p">,</span><span class="n">n</span><span class="p">])</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">minmn</span><span class="p">):</span>
            <span class="k">if</span> <span class="p">(</span><span class="n">pivot</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">):</span>
                <span class="c">## Bring the column of largest norm into the pivot position</span>
                <span class="n">rmax</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">max</span><span class="p">(</span><span class="n">rdiag</span><span class="p">[</span><span class="n">j</span><span class="p">:])</span>
                <span class="n">kmax</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">(</span><span class="n">rdiag</span><span class="p">[</span><span class="n">j</span><span class="p">:]</span> <span class="o">==</span> <span class="n">rmax</span><span class="p">))[</span><span class="mi">0</span><span class="p">]</span>
                <span class="n">ct</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">kmax</span><span class="p">)</span>
                <span class="n">kmax</span> <span class="o">=</span> <span class="n">kmax</span> <span class="o">+</span> <span class="n">j</span>
                <span class="k">if</span> <span class="n">ct</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
                    <span class="n">kmax</span> <span class="o">=</span> <span class="n">kmax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>

                    <span class="c">## Exchange rows via the pivot only.  Avoid actually exchanging</span>
                    <span class="c">## the rows, in case there is lots of memory transfer.  The</span>
                    <span class="c">## exchange occurs later, within the body of MPFIT, after the</span>
                    <span class="c">## extraneous columns of the matrix have been shed.</span>
                    <span class="k">if</span> <span class="n">kmax</span> <span class="o">!=</span> <span class="n">j</span><span class="p">:</span>
                        <span class="n">temp</span> <span class="o">=</span> <span class="n">ipvt</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="p">;</span> <span class="n">ipvt</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">ipvt</span><span class="p">[</span><span class="n">kmax</span><span class="p">]</span> <span class="p">;</span> <span class="n">ipvt</span><span class="p">[</span><span class="n">kmax</span><span class="p">]</span> <span class="o">=</span> <span class="n">temp</span>
                        <span class="n">rdiag</span><span class="p">[</span><span class="n">kmax</span><span class="p">]</span> <span class="o">=</span> <span class="n">rdiag</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
                        <span class="n">wa</span><span class="p">[</span><span class="n">kmax</span><span class="p">]</span> <span class="o">=</span> <span class="n">wa</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>

            <span class="c">## Compute the householder transformation to reduce the jth</span>
            <span class="c">## column of A to a multiple of the jth unit vector</span>
            <span class="n">lj</span> <span class="o">=</span> <span class="n">ipvt</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
            <span class="n">ajj</span> <span class="o">=</span> <span class="n">a</span><span class="p">[</span><span class="n">j</span><span class="p">:,</span><span class="n">lj</span><span class="p">]</span>
            <span class="n">ajnorm</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">enorm</span><span class="p">(</span><span class="n">ajj</span><span class="p">)</span>
            <span class="k">if</span> <span class="n">ajnorm</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span> <span class="k">break</span>
            <span class="k">if</span> <span class="n">a</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">&lt;</span> <span class="mi">0</span><span class="p">:</span> <span class="n">ajnorm</span> <span class="o">=</span> <span class="o">-</span><span class="n">ajnorm</span>

            <span class="n">ajj</span> <span class="o">=</span> <span class="n">ajj</span> <span class="o">/</span> <span class="n">ajnorm</span>
            <span class="n">ajj</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">ajj</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">+</span> <span class="mi">1</span>
            <span class="c">## *** Note optimization a(j:*,j)</span>
            <span class="n">a</span><span class="p">[</span><span class="n">j</span><span class="p">:,</span><span class="n">lj</span><span class="p">]</span> <span class="o">=</span> <span class="n">ajj</span>

            <span class="c">## Apply the transformation to the remaining columns</span>
            <span class="c">## and update the norms</span>

            <span class="c">## NOTE to SELF: tried to optimize this by removing the loop,</span>
            <span class="c">## but it actually got slower.  Reverted to &quot;for&quot; loop to keep</span>
            <span class="c">## it simple.</span>
            <span class="k">if</span> <span class="p">(</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span> <span class="o">&lt;</span> <span class="n">n</span><span class="p">):</span>
                <span class="k">for</span> <span class="n">k</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span> <span class="n">n</span><span class="p">):</span>
                    <span class="n">lk</span> <span class="o">=</span> <span class="n">ipvt</span><span class="p">[</span><span class="n">k</span><span class="p">]</span>
                    <span class="n">ajk</span> <span class="o">=</span> <span class="n">a</span><span class="p">[</span><span class="n">j</span><span class="p">:,</span><span class="n">lk</span><span class="p">]</span>
                    <span class="c">## *** Note optimization a(j:*,lk)</span>
                    <span class="c">## (corrected 20 Jul 2000)</span>
                    <span class="k">if</span> <span class="n">a</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">lj</span><span class="p">]</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">:</span>
                        <span class="n">a</span><span class="p">[</span><span class="n">j</span><span class="p">:,</span><span class="n">lk</span><span class="p">]</span> <span class="o">=</span> <span class="n">ajk</span> <span class="o">-</span> <span class="n">ajj</span> <span class="o">*</span> <span class="nb">sum</span><span class="p">(</span><span class="n">ajk</span><span class="o">*</span><span class="n">ajj</span><span class="p">)</span><span class="o">/</span><span class="n">a</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">lj</span><span class="p">]</span>
                        <span class="k">if</span> <span class="p">((</span><span class="n">pivot</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">rdiag</span><span class="p">[</span><span class="n">k</span><span class="p">]</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">)):</span>
                            <span class="n">temp</span> <span class="o">=</span> <span class="n">a</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">lk</span><span class="p">]</span><span class="o">/</span><span class="n">rdiag</span><span class="p">[</span><span class="n">k</span><span class="p">]</span>
                            <span class="n">rdiag</span><span class="p">[</span><span class="n">k</span><span class="p">]</span> <span class="o">=</span> <span class="n">rdiag</span><span class="p">[</span><span class="n">k</span><span class="p">]</span> <span class="o">*</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">numpy</span><span class="o">.</span><span class="n">max</span><span class="p">([(</span><span class="mf">1.</span><span class="o">-</span><span class="n">temp</span><span class="o">**</span><span class="mi">2</span><span class="p">),</span> <span class="mf">0.</span><span class="p">]))</span>
                            <span class="n">temp</span> <span class="o">=</span> <span class="n">rdiag</span><span class="p">[</span><span class="n">k</span><span class="p">]</span><span class="o">/</span><span class="n">wa</span><span class="p">[</span><span class="n">k</span><span class="p">]</span>
                            <span class="k">if</span> <span class="p">((</span><span class="mf">0.05</span><span class="o">*</span><span class="n">temp</span><span class="o">*</span><span class="n">temp</span><span class="p">)</span> <span class="o">&lt;=</span> <span class="n">machep</span><span class="p">):</span>
                                <span class="n">rdiag</span><span class="p">[</span><span class="n">k</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">enorm</span><span class="p">(</span><span class="n">a</span><span class="p">[</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">:,</span><span class="n">lk</span><span class="p">])</span>
                                <span class="n">wa</span><span class="p">[</span><span class="n">k</span><span class="p">]</span> <span class="o">=</span> <span class="n">rdiag</span><span class="p">[</span><span class="n">k</span><span class="p">]</span>
            <span class="n">rdiag</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="o">-</span><span class="n">ajnorm</span>
        <span class="k">return</span><span class="p">([</span><span class="n">a</span><span class="p">,</span> <span class="n">ipvt</span><span class="p">,</span> <span class="n">rdiag</span><span class="p">,</span> <span class="n">acnorm</span><span class="p">])</span>

    
    <span class="c">#    Original FORTRAN documentation</span>
    <span class="c">#    **********</span>
    <span class="c">#</span>
    <span class="c">#    subroutine qrsolv</span>
    <span class="c">#</span>
    <span class="c">#    given an m by n matrix a, an n by n diagonal matrix d,</span>
    <span class="c">#    and an m-vector b, the problem is to determine an x which</span>
    <span class="c">#    solves the system</span>
    <span class="c">#</span>
    <span class="c">#          a*x = b ,     d*x = 0 ,</span>
    <span class="c">#</span>
    <span class="c">#    in the least squares sense.</span>
    <span class="c">#</span>
    <span class="c">#    this subroutine completes the solution of the problem</span>
    <span class="c">#    if it is provided with the necessary information from the</span>
    <span class="c">#    factorization, with column pivoting, of a. that is, if</span>
    <span class="c">#    a*p = q*r, where p is a permutation matrix, q has orthogonal</span>
    <span class="c">#    columns, and r is an upper triangular matrix with diagonal</span>
    <span class="c">#    elements of nonincreasing magnitude, then qrsolv expects</span>
    <span class="c">#    the full upper triangle of r, the permutation matrix p,</span>
    <span class="c">#    and the first n components of (q transpose)*b. the system</span>
    <span class="c">#    a*x = b, d*x = 0, is then equivalent to</span>
    <span class="c">#</span>
    <span class="c">#                 t    t</span>
    <span class="c">#          r*z = q *b ,  p *d*p*z = 0 ,</span>
    <span class="c">#</span>
    <span class="c">#    where x = p*z. if this system does not have full rank,</span>
    <span class="c">#    then a least squares solution is obtained. on output qrsolv</span>
    <span class="c">#    also provides an upper triangular matrix s such that</span>
    <span class="c">#</span>
    <span class="c">#           t   t              t</span>
    <span class="c">#          p *(a *a + d*d)*p = s *s .</span>
    <span class="c">#</span>
    <span class="c">#    s is computed within qrsolv and may be of separate interest.</span>
    <span class="c">#</span>
    <span class="c">#    the subroutine statement is</span>
    <span class="c">#</span>
    <span class="c">#      subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa)</span>
    <span class="c">#</span>
    <span class="c">#    where</span>
    <span class="c">#</span>
    <span class="c">#      n is a positive integer input variable set to the order of r.</span>
    <span class="c">#</span>
    <span class="c">#      r is an n by n array. on input the full upper triangle</span>
    <span class="c">#        must contain the full upper triangle of the matrix r.</span>
    <span class="c">#        on output the full upper triangle is unaltered, and the</span>
    <span class="c">#        strict lower triangle contains the strict upper triangle</span>
    <span class="c">#        (transposed) of the upper triangular matrix s.</span>
    <span class="c">#</span>
    <span class="c">#      ldr is a positive integer input variable not less than n</span>
    <span class="c">#        which specifies the leading dimension of the array r.</span>
    <span class="c">#</span>
    <span class="c">#      ipvt is an integer input array of length n which defines the</span>
    <span class="c">#        permutation matrix p such that a*p = q*r. column j of p</span>
    <span class="c">#        is column ipvt(j) of the identity matrix.</span>
    <span class="c">#</span>
    <span class="c">#      diag is an input array of length n which must contain the</span>
    <span class="c">#        diagonal elements of the matrix d.</span>
    <span class="c">#</span>
    <span class="c">#      qtb is an input array of length n which must contain the first</span>
    <span class="c">#        n elements of the vector (q transpose)*b.</span>
    <span class="c">#</span>
    <span class="c">#      x is an output array of length n which contains the least</span>
    <span class="c">#        squares solution of the system a*x = b, d*x = 0.</span>
    <span class="c">#</span>
    <span class="c">#      sdiag is an output array of length n which contains the</span>
    <span class="c">#        diagonal elements of the upper triangular matrix s.</span>
    <span class="c">#</span>
    <span class="c">#      wa is a work array of length n.</span>
    <span class="c">#</span>
    <span class="c">#    subprograms called</span>
    <span class="c">#</span>
    <span class="c">#      fortran-supplied ... dabs,dsqrt</span>
    <span class="c">#</span>
    <span class="c">#    argonne national laboratory. minpack project. march 1980.</span>
    <span class="c">#    burton s. garbow, kenneth e. hillstrom, jorge j. more</span>
    <span class="c">#</span>
    </div>
<div class="viewcode-block" id="mpfit.qrsolv"><a class="viewcode-back" href="../../agpy.html#agpy.mpfit_custom.mpfit.qrsolv">[docs]</a>    <span class="k">def</span> <span class="nf">qrsolv</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">r</span><span class="p">,</span> <span class="n">ipvt</span><span class="p">,</span> <span class="n">diag</span><span class="p">,</span> <span class="n">qtb</span><span class="p">,</span> <span class="n">sdiag</span><span class="p">):</span>
        <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">debug</span><span class="p">):</span> <span class="k">print</span> <span class="s">&#39;Entering qrsolv...&#39;</span>
        <span class="n">sz</span> <span class="o">=</span> <span class="n">shape</span><span class="p">(</span><span class="n">r</span><span class="p">)</span>
        <span class="n">m</span> <span class="o">=</span> <span class="n">sz</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
        <span class="n">n</span> <span class="o">=</span> <span class="n">sz</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>

        <span class="c">## copy r and (q transpose)*b to preserve input and initialize s.</span>
        <span class="c">## in particular, save the diagonal elements of r in x.</span>

        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
            <span class="n">r</span><span class="p">[</span><span class="n">j</span><span class="p">:</span><span class="n">n</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">j</span><span class="p">:</span><span class="n">n</span><span class="p">]</span>
        <span class="n">x</span> <span class="o">=</span> <span class="n">diagonal</span><span class="p">(</span><span class="n">r</span><span class="p">)</span>
        <span class="n">wa</span> <span class="o">=</span> <span class="n">qtb</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>

        <span class="c">## Eliminate the diagonal matrix d using a givens rotation</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
            <span class="n">l</span> <span class="o">=</span> <span class="n">ipvt</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
            <span class="k">if</span> <span class="p">(</span><span class="n">diag</span><span class="p">[</span><span class="n">l</span><span class="p">]</span> <span class="o">==</span> <span class="mi">0</span><span class="p">):</span> <span class="k">break</span>
            <span class="n">sdiag</span><span class="p">[</span><span class="n">j</span><span class="p">:]</span> <span class="o">=</span> <span class="mi">0</span>
            <span class="n">sdiag</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">diag</span><span class="p">[</span><span class="n">l</span><span class="p">]</span>

            <span class="c">## The transformations to eliminate the row of d modify only a</span>
            <span class="c">## single element of (q transpose)*b beyond the first n, which</span>
            <span class="c">## is initially zero.</span>

            <span class="n">qtbpj</span> <span class="o">=</span> <span class="mf">0.</span>
            <span class="k">for</span> <span class="n">k</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">j</span><span class="p">,</span><span class="n">n</span><span class="p">):</span>
                <span class="k">if</span> <span class="p">(</span><span class="n">sdiag</span><span class="p">[</span><span class="n">k</span><span class="p">]</span> <span class="o">==</span> <span class="mi">0</span><span class="p">):</span> <span class="k">break</span>
                <span class="k">if</span> <span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="n">k</span><span class="p">,</span><span class="n">k</span><span class="p">])</span> <span class="o">&lt;</span> <span class="nb">abs</span><span class="p">(</span><span class="n">sdiag</span><span class="p">[</span><span class="n">k</span><span class="p">])):</span>
                    <span class="n">cotan</span>  <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="n">k</span><span class="p">,</span><span class="n">k</span><span class="p">]</span><span class="o">/</span><span class="n">sdiag</span><span class="p">[</span><span class="n">k</span><span class="p">]</span>
                    <span class="n">sine</span>   <span class="o">=</span> <span class="mf">0.5</span><span class="o">/</span><span class="n">sqrt</span><span class="p">(</span><span class="o">.</span><span class="mi">25</span> <span class="o">+</span> <span class="o">.</span><span class="mi">25</span><span class="o">*</span><span class="n">cotan</span><span class="o">*</span><span class="n">cotan</span><span class="p">)</span>
                    <span class="n">cosine</span> <span class="o">=</span> <span class="n">sine</span><span class="o">*</span><span class="n">cotan</span>
                <span class="k">else</span><span class="p">:</span>
                    <span class="n">tang</span>   <span class="o">=</span> <span class="n">sdiag</span><span class="p">[</span><span class="n">k</span><span class="p">]</span><span class="o">/</span><span class="n">r</span><span class="p">[</span><span class="n">k</span><span class="p">,</span><span class="n">k</span><span class="p">]</span>
                    <span class="n">cosine</span> <span class="o">=</span> <span class="mf">0.5</span><span class="o">/</span><span class="n">sqrt</span><span class="p">(</span><span class="o">.</span><span class="mi">25</span> <span class="o">+</span> <span class="o">.</span><span class="mi">25</span><span class="o">*</span><span class="n">tang</span><span class="o">*</span><span class="n">tang</span><span class="p">)</span>
                    <span class="n">sine</span>   <span class="o">=</span> <span class="n">cosine</span><span class="o">*</span><span class="n">tang</span>

                <span class="c">## Compute the modified diagonal element of r and the</span>
                <span class="c">## modified element of ((q transpose)*b,0).</span>
                <span class="n">r</span><span class="p">[</span><span class="n">k</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">=</span> <span class="n">cosine</span><span class="o">*</span><span class="n">r</span><span class="p">[</span><span class="n">k</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">+</span> <span class="n">sine</span><span class="o">*</span><span class="n">sdiag</span><span class="p">[</span><span class="n">k</span><span class="p">]</span>
                <span class="n">temp</span> <span class="o">=</span> <span class="n">cosine</span><span class="o">*</span><span class="n">wa</span><span class="p">[</span><span class="n">k</span><span class="p">]</span> <span class="o">+</span> <span class="n">sine</span><span class="o">*</span><span class="n">qtbpj</span>
                <span class="n">qtbpj</span> <span class="o">=</span> <span class="o">-</span><span class="n">sine</span><span class="o">*</span><span class="n">wa</span><span class="p">[</span><span class="n">k</span><span class="p">]</span> <span class="o">+</span> <span class="n">cosine</span><span class="o">*</span><span class="n">qtbpj</span>
                <span class="n">wa</span><span class="p">[</span><span class="n">k</span><span class="p">]</span> <span class="o">=</span> <span class="n">temp</span>

                <span class="c">## Accumulate the transformation in the row of s</span>
                <span class="k">if</span> <span class="p">(</span><span class="n">n</span> <span class="o">&gt;</span> <span class="n">k</span><span class="o">+</span><span class="mi">1</span><span class="p">):</span>
                    <span class="n">temp</span> <span class="o">=</span> <span class="n">cosine</span><span class="o">*</span><span class="n">r</span><span class="p">[</span><span class="n">k</span><span class="o">+</span><span class="mi">1</span><span class="p">:</span><span class="n">n</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">+</span> <span class="n">sine</span><span class="o">*</span><span class="n">sdiag</span><span class="p">[</span><span class="n">k</span><span class="o">+</span><span class="mi">1</span><span class="p">:</span><span class="n">n</span><span class="p">]</span>
                    <span class="n">sdiag</span><span class="p">[</span><span class="n">k</span><span class="o">+</span><span class="mi">1</span><span class="p">:</span><span class="n">n</span><span class="p">]</span> <span class="o">=</span> <span class="o">-</span><span class="n">sine</span><span class="o">*</span><span class="n">r</span><span class="p">[</span><span class="n">k</span><span class="o">+</span><span class="mi">1</span><span class="p">:</span><span class="n">n</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">+</span> <span class="n">cosine</span><span class="o">*</span><span class="n">sdiag</span><span class="p">[</span><span class="n">k</span><span class="o">+</span><span class="mi">1</span><span class="p">:</span><span class="n">n</span><span class="p">]</span>
                    <span class="n">r</span><span class="p">[</span><span class="n">k</span><span class="o">+</span><span class="mi">1</span><span class="p">:</span><span class="n">n</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">=</span> <span class="n">temp</span>
            <span class="n">sdiag</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">j</span><span class="p">]</span>
            <span class="n">r</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>

        <span class="c">## Solve the triangular system for z.  If the system is singular</span>
        <span class="c">## then obtain a least squares solution</span>
        <span class="n">nsing</span> <span class="o">=</span> <span class="n">n</span>
        <span class="n">wh</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">(</span><span class="n">sdiag</span> <span class="o">==</span> <span class="mi">0</span><span class="p">))[</span><span class="mi">0</span><span class="p">]</span>
        <span class="k">if</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">wh</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">):</span>
            <span class="n">nsing</span> <span class="o">=</span> <span class="n">wh</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
            <span class="n">wa</span><span class="p">[</span><span class="n">nsing</span><span class="p">:]</span> <span class="o">=</span> <span class="mi">0</span>

        <span class="k">if</span> <span class="p">(</span><span class="n">nsing</span> <span class="o">&gt;=</span> <span class="mi">1</span><span class="p">):</span>
            <span class="n">wa</span><span class="p">[</span><span class="n">nsing</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">wa</span><span class="p">[</span><span class="n">nsing</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">/</span><span class="n">sdiag</span><span class="p">[</span><span class="n">nsing</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="c">## Degenerate case</span>
            <span class="c">## *** Reverse loop ***</span>
            <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nsing</span><span class="o">-</span><span class="mi">2</span><span class="p">,</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
                <span class="n">sum0</span> <span class="o">=</span> <span class="nb">sum</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">:</span><span class="n">nsing</span><span class="p">,</span><span class="n">j</span><span class="p">]</span><span class="o">*</span><span class="n">wa</span><span class="p">[</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">:</span><span class="n">nsing</span><span class="p">])</span>
                <span class="n">wa</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">wa</span><span class="p">[</span><span class="n">j</span><span class="p">]</span><span class="o">-</span><span class="n">sum0</span><span class="p">)</span><span class="o">/</span><span class="n">sdiag</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>

        <span class="c">## Permute the components of z back to components of x</span>
        <span class="n">put</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">ipvt</span><span class="p">,</span> <span class="n">wa</span><span class="p">)</span>
        <span class="k">return</span><span class="p">(</span><span class="n">r</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">sdiag</span><span class="p">)</span>



    
    <span class="c">#    Original FORTRAN documentation</span>
    <span class="c">#</span>
    <span class="c">#    subroutine lmpar</span>
    <span class="c">#</span>
    <span class="c">#    given an m by n matrix a, an n by n nonsingular diagonal</span>
    <span class="c">#    matrix d, an m-vector b, and a positive number delta,</span>
    <span class="c">#    the problem is to determine a value for the parameter</span>
    <span class="c">#    par such that if x solves the system</span>
    <span class="c">#</span>
    <span class="c">#       a*x = b ,    sqrt(par)*d*x = 0 ,</span>
    <span class="c">#</span>
    <span class="c">#    in the least squares sense, and dxnorm is the euclidean</span>
    <span class="c">#    norm of d*x, then either par is zero and</span>
    <span class="c">#</span>
    <span class="c">#       (dxnorm-delta) .le. 0.1*delta ,</span>
    <span class="c">#</span>
    <span class="c">#    or par is positive and</span>
    <span class="c">#</span>
    <span class="c">#       abs(dxnorm-delta) .le. 0.1*delta .</span>
    <span class="c">#</span>
    <span class="c">#    this subroutine completes the solution of the problem</span>
    <span class="c">#    if it is provided with the necessary information from the</span>
    <span class="c">#    qr factorization, with column pivoting, of a. that is, if</span>
    <span class="c">#    a*p = q*r, where p is a permutation matrix, q has orthogonal</span>
    <span class="c">#    columns, and r is an upper triangular matrix with diagonal</span>
    <span class="c">#    elements of nonincreasing magnitude, then lmpar expects</span>
    <span class="c">#    the full upper triangle of r, the permutation matrix p,</span>
    <span class="c">#    and the first n components of (q transpose)*b. on output</span>
    <span class="c">#    lmpar also provides an upper triangular matrix s such that</span>
    <span class="c">#</span>
    <span class="c">#        t   t                 t</span>
    <span class="c">#       p *(a *a + par*d*d)*p = s *s .</span>
    <span class="c">#</span>
    <span class="c">#    s is employed within lmpar and may be of separate interest.</span>
    <span class="c">#</span>
    <span class="c">#    only a few iterations are generally needed for convergence</span>
    <span class="c">#    of the algorithm. if, however, the limit of 10 iterations</span>
    <span class="c">#    is reached, then the output par will contain the best</span>
    <span class="c">#    value obtained so far.</span>
    <span class="c">#</span>
    <span class="c">#    the subroutine statement is</span>
    <span class="c">#</span>
    <span class="c">#   subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag,</span>
    <span class="c">#                    wa1,wa2)</span>
    <span class="c">#</span>
    <span class="c">#    where</span>
    <span class="c">#</span>
    <span class="c">#   n is a positive integer input variable set to the order of r.</span>
    <span class="c">#</span>
    <span class="c">#   r is an n by n array. on input the full upper triangle</span>
    <span class="c">#     must contain the full upper triangle of the matrix r.</span>
    <span class="c">#     on output the full upper triangle is unaltered, and the</span>
    <span class="c">#     strict lower triangle contains the strict upper triangle</span>
    <span class="c">#     (transposed) of the upper triangular matrix s.</span>
    <span class="c">#</span>
    <span class="c">#   ldr is a positive integer input variable not less than n</span>
    <span class="c">#     which specifies the leading dimension of the array r.</span>
    <span class="c">#</span>
    <span class="c">#   ipvt is an integer input array of length n which defines the</span>
    <span class="c">#     permutation matrix p such that a*p = q*r. column j of p</span>
    <span class="c">#     is column ipvt(j) of the identity matrix.</span>
    <span class="c">#</span>
    <span class="c">#   diag is an input array of length n which must contain the</span>
    <span class="c">#     diagonal elements of the matrix d.</span>
    <span class="c">#</span>
    <span class="c">#   qtb is an input array of length n which must contain the first</span>
    <span class="c">#     n elements of the vector (q transpose)*b.</span>
    <span class="c">#</span>
    <span class="c">#   delta is a positive input variable which specifies an upper</span>
    <span class="c">#     bound on the euclidean norm of d*x.</span>
    <span class="c">#</span>
    <span class="c">#   par is a nonnegative variable. on input par contains an</span>
    <span class="c">#     initial estimate of the levenberg-marquardt parameter.</span>
    <span class="c">#     on output par contains the final estimate.</span>
    <span class="c">#</span>
    <span class="c">#   x is an output array of length n which contains the least</span>
    <span class="c">#     squares solution of the system a*x = b, sqrt(par)*d*x = 0,</span>
    <span class="c">#     for the output par.</span>
    <span class="c">#</span>
    <span class="c">#   sdiag is an output array of length n which contains the</span>
    <span class="c">#     diagonal elements of the upper triangular matrix s.</span>
    <span class="c">#</span>
    <span class="c">#   wa1 and wa2 are work arrays of length n.</span>
    <span class="c">#</span>
    <span class="c">#    subprograms called</span>
    <span class="c">#</span>
    <span class="c">#   minpack-supplied ... dpmpar,enorm,qrsolv</span>
    <span class="c">#</span>
    <span class="c">#   fortran-supplied ... dabs,dmax1,dmin1,dsqrt</span>
    <span class="c">#</span>
    <span class="c">#    argonne national laboratory. minpack project. march 1980.</span>
    <span class="c">#    burton s. garbow, kenneth e. hillstrom, jorge j. more</span>
    <span class="c">#</span>
    </div>
<div class="viewcode-block" id="mpfit.lmpar"><a class="viewcode-back" href="../../agpy.html#agpy.mpfit_custom.mpfit.lmpar">[docs]</a>    <span class="k">def</span> <span class="nf">lmpar</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">r</span><span class="p">,</span> <span class="n">ipvt</span><span class="p">,</span> <span class="n">diag</span><span class="p">,</span> <span class="n">qtb</span><span class="p">,</span> <span class="n">delta</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">sdiag</span><span class="p">,</span> <span class="n">par</span><span class="o">=</span><span class="bp">None</span><span class="p">):</span>

        <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">debug</span><span class="p">):</span> <span class="k">print</span> <span class="s">&#39;Entering lmpar...&#39;</span>
        <span class="n">dwarf</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">machar</span><span class="o">.</span><span class="n">minnum</span>
        <span class="n">sz</span> <span class="o">=</span> <span class="n">shape</span><span class="p">(</span><span class="n">r</span><span class="p">)</span>
        <span class="n">m</span> <span class="o">=</span> <span class="n">sz</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
        <span class="n">n</span> <span class="o">=</span> <span class="n">sz</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>

        <span class="c">## Compute and store in x the gauss-newton direction.  If the</span>
        <span class="c">## jacobian is rank-deficient, obtain a least-squares solution</span>
        <span class="n">nsing</span> <span class="o">=</span> <span class="n">n</span>
        <span class="n">wa1</span> <span class="o">=</span> <span class="n">qtb</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
        <span class="n">wh</span> <span class="o">=</span> <span class="p">(</span><span class="n">nonzero</span><span class="p">(</span><span class="n">diagonal</span><span class="p">(</span><span class="n">r</span><span class="p">)</span> <span class="o">==</span> <span class="mi">0</span><span class="p">))[</span><span class="mi">0</span><span class="p">]</span>
        <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">wh</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
            <span class="n">nsing</span> <span class="o">=</span> <span class="n">wh</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
            <span class="n">wa1</span><span class="p">[</span><span class="n">wh</span><span class="p">[</span><span class="mi">0</span><span class="p">]:]</span> <span class="o">=</span> <span class="mi">0</span>
        <span class="k">if</span> <span class="n">nsing</span> <span class="o">&gt;=</span> <span class="mi">1</span><span class="p">:</span>
            <span class="c">## *** Reverse loop ***</span>
            <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nsing</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
                <span class="n">wa1</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">wa1</span><span class="p">[</span><span class="n">j</span><span class="p">]</span><span class="o">/</span><span class="n">r</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">j</span><span class="p">]</span>
                <span class="k">if</span> <span class="p">(</span><span class="n">j</span><span class="o">-</span><span class="mi">1</span> <span class="o">&gt;=</span> <span class="mi">0</span><span class="p">):</span>
                    <span class="n">wa1</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">wa1</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="p">]</span> <span class="o">-</span> <span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="p">,</span><span class="n">j</span><span class="p">]</span><span class="o">*</span><span class="n">wa1</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>

        <span class="c">## Note: ipvt here is a permutation array</span>
        <span class="n">put</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">ipvt</span><span class="p">,</span> <span class="n">wa1</span><span class="p">)</span>

        <span class="c">## Initialize the iteration counter.  Evaluate the function at the</span>
        <span class="c">## origin, and test for acceptance of the gauss-newton direction</span>
        <span class="nb">iter</span> <span class="o">=</span> <span class="mi">0</span>
        <span class="n">wa2</span> <span class="o">=</span> <span class="n">diag</span> <span class="o">*</span> <span class="n">x</span>
        <span class="n">dxnorm</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">enorm</span><span class="p">(</span><span class="n">wa2</span><span class="p">)</span>
        <span class="n">fp</span> <span class="o">=</span> <span class="n">dxnorm</span> <span class="o">-</span> <span class="n">delta</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">fp</span> <span class="o">&lt;=</span> <span class="mf">0.1</span><span class="o">*</span><span class="n">delta</span><span class="p">):</span>
            <span class="k">return</span><span class="p">[</span><span class="n">r</span><span class="p">,</span> <span class="mf">0.</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">sdiag</span><span class="p">]</span>

        <span class="c">## If the jacobian is not rank deficient, the newton step provides a</span>
        <span class="c">## lower bound, parl, for the zero of the function.  Otherwise set</span>
        <span class="c">## this bound to zero.</span>

        <span class="n">parl</span> <span class="o">=</span> <span class="mf">0.</span>
        <span class="k">if</span> <span class="n">nsing</span> <span class="o">&gt;=</span> <span class="n">n</span><span class="p">:</span>
            <span class="n">wa1</span> <span class="o">=</span> <span class="n">take</span><span class="p">(</span><span class="n">diag</span><span class="p">,</span> <span class="n">ipvt</span><span class="p">)</span><span class="o">*</span><span class="n">take</span><span class="p">(</span><span class="n">wa2</span><span class="p">,</span> <span class="n">ipvt</span><span class="p">)</span><span class="o">/</span><span class="n">dxnorm</span>
            <span class="n">wa1</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">wa1</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">/</span> <span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span> <span class="c">## Degenerate case</span>
            <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span><span class="n">n</span><span class="p">):</span>   <span class="c">## Note &quot;1&quot; here, not zero</span>
                <span class="n">sum0</span> <span class="o">=</span> <span class="nb">sum</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="p">,</span><span class="n">j</span><span class="p">]</span><span class="o">*</span><span class="n">wa1</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="p">])</span>
                <span class="n">wa1</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">wa1</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">-</span> <span class="n">sum0</span><span class="p">)</span><span class="o">/</span><span class="n">r</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">j</span><span class="p">]</span>

            <span class="n">temp</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">enorm</span><span class="p">(</span><span class="n">wa1</span><span class="p">)</span>
            <span class="n">parl</span> <span class="o">=</span> <span class="p">((</span><span class="n">fp</span><span class="o">/</span><span class="n">delta</span><span class="p">)</span><span class="o">/</span><span class="n">temp</span><span class="p">)</span><span class="o">/</span><span class="n">temp</span>

        <span class="c">## Calculate an upper bound, paru, for the zero of the function</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
            <span class="n">sum0</span> <span class="o">=</span> <span class="nb">sum</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span><span class="n">j</span><span class="p">]</span><span class="o">*</span><span class="n">qtb</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">])</span>
            <span class="n">wa1</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">sum0</span><span class="o">/</span><span class="n">diag</span><span class="p">[</span><span class="n">ipvt</span><span class="p">[</span><span class="n">j</span><span class="p">]]</span>
        <span class="n">gnorm</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">enorm</span><span class="p">(</span><span class="n">wa1</span><span class="p">)</span>
        <span class="n">paru</span> <span class="o">=</span> <span class="n">gnorm</span><span class="o">/</span><span class="n">delta</span>
        <span class="k">if</span> <span class="n">paru</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span> <span class="n">paru</span> <span class="o">=</span> <span class="n">dwarf</span><span class="o">/</span><span class="n">numpy</span><span class="o">.</span><span class="n">min</span><span class="p">([</span><span class="n">delta</span><span class="p">,</span><span class="mf">0.1</span><span class="p">])</span>

        <span class="c">## If the input par lies outside of the interval (parl,paru), set</span>
        <span class="c">## par to the closer endpoint</span>

        <span class="n">par</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">max</span><span class="p">([</span><span class="n">par</span><span class="p">,</span><span class="n">parl</span><span class="p">])</span>
        <span class="n">par</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">min</span><span class="p">([</span><span class="n">par</span><span class="p">,</span><span class="n">paru</span><span class="p">])</span>
        <span class="k">if</span> <span class="n">par</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span> <span class="n">par</span> <span class="o">=</span> <span class="n">gnorm</span><span class="o">/</span><span class="n">dxnorm</span>

        <span class="c">## Beginning of an interation</span>
        <span class="k">while</span><span class="p">(</span><span class="mi">1</span><span class="p">):</span>
            <span class="nb">iter</span> <span class="o">=</span> <span class="nb">iter</span> <span class="o">+</span> <span class="mi">1</span>

            <span class="c">## Evaluate the function at the current value of par</span>
            <span class="k">if</span> <span class="n">par</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span> <span class="n">par</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">max</span><span class="p">([</span><span class="n">dwarf</span><span class="p">,</span> <span class="n">paru</span><span class="o">*</span><span class="mf">0.001</span><span class="p">])</span>
            <span class="n">temp</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">par</span><span class="p">)</span>
            <span class="n">wa1</span> <span class="o">=</span> <span class="n">temp</span> <span class="o">*</span> <span class="n">diag</span>
            <span class="p">[</span><span class="n">r</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">sdiag</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">qrsolv</span><span class="p">(</span><span class="n">r</span><span class="p">,</span> <span class="n">ipvt</span><span class="p">,</span> <span class="n">wa1</span><span class="p">,</span> <span class="n">qtb</span><span class="p">,</span> <span class="n">sdiag</span><span class="p">)</span>
            <span class="n">wa2</span> <span class="o">=</span> <span class="n">diag</span><span class="o">*</span><span class="n">x</span>
            <span class="n">dxnorm</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">enorm</span><span class="p">(</span><span class="n">wa2</span><span class="p">)</span>
            <span class="n">temp</span> <span class="o">=</span> <span class="n">fp</span>
            <span class="n">fp</span> <span class="o">=</span> <span class="n">dxnorm</span> <span class="o">-</span> <span class="n">delta</span>

            <span class="k">if</span> <span class="p">((</span><span class="nb">abs</span><span class="p">(</span><span class="n">fp</span><span class="p">)</span> <span class="o">&lt;=</span> <span class="mf">0.1</span><span class="o">*</span><span class="n">delta</span><span class="p">)</span> <span class="ow">or</span>
               <span class="p">((</span><span class="n">parl</span> <span class="o">==</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">fp</span> <span class="o">&lt;=</span> <span class="n">temp</span><span class="p">)</span> <span class="ow">and</span> <span class="p">(</span><span class="n">temp</span> <span class="o">&lt;</span> <span class="mi">0</span><span class="p">))</span> <span class="ow">or</span>
               <span class="p">(</span><span class="nb">iter</span> <span class="o">==</span> <span class="mi">10</span><span class="p">)):</span> <span class="k">break</span><span class="p">;</span>

            <span class="c">## Compute the newton correction</span>
            <span class="n">wa1</span> <span class="o">=</span> <span class="n">take</span><span class="p">(</span><span class="n">diag</span><span class="p">,</span> <span class="n">ipvt</span><span class="p">)</span><span class="o">*</span><span class="n">take</span><span class="p">(</span><span class="n">wa2</span><span class="p">,</span> <span class="n">ipvt</span><span class="p">)</span><span class="o">/</span><span class="n">dxnorm</span>

            <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
                <span class="n">wa1</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">wa1</span><span class="p">[</span><span class="n">j</span><span class="p">]</span><span class="o">/</span><span class="n">sdiag</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
                <span class="n">wa1</span><span class="p">[</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">:</span><span class="n">n</span><span class="p">]</span> <span class="o">=</span> <span class="n">wa1</span><span class="p">[</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">:</span><span class="n">n</span><span class="p">]</span> <span class="o">-</span> <span class="n">r</span><span class="p">[</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">:</span><span class="n">n</span><span class="p">,</span><span class="n">j</span><span class="p">]</span><span class="o">*</span><span class="n">wa1</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
            <span class="n">wa1</span><span class="p">[</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">wa1</span><span class="p">[</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">/</span><span class="n">sdiag</span><span class="p">[</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="c">## Degenerate case</span>

            <span class="n">temp</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">enorm</span><span class="p">(</span><span class="n">wa1</span><span class="p">)</span>
            <span class="n">parc</span> <span class="o">=</span> <span class="p">((</span><span class="n">fp</span><span class="o">/</span><span class="n">delta</span><span class="p">)</span><span class="o">/</span><span class="n">temp</span><span class="p">)</span><span class="o">/</span><span class="n">temp</span>

            <span class="c">## Depending on the sign of the function, update parl or paru</span>
            <span class="k">if</span> <span class="n">fp</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span> <span class="n">parl</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">max</span><span class="p">([</span><span class="n">parl</span><span class="p">,</span><span class="n">par</span><span class="p">])</span>
            <span class="k">if</span> <span class="n">fp</span> <span class="o">&lt;</span> <span class="mi">0</span><span class="p">:</span> <span class="n">paru</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">min</span><span class="p">([</span><span class="n">paru</span><span class="p">,</span><span class="n">par</span><span class="p">])</span>

            <span class="c">## Compute an improved estimate for par</span>
            <span class="n">par</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">max</span><span class="p">([</span><span class="n">parl</span><span class="p">,</span> <span class="n">par</span><span class="o">+</span><span class="n">parc</span><span class="p">])</span>

            <span class="c">## End of an iteration</span>

        <span class="c">## Termination</span>
        <span class="k">return</span><span class="p">[</span><span class="n">r</span><span class="p">,</span> <span class="n">par</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">sdiag</span><span class="p">]</span>

    
    <span class="c">## Procedure to tie one parameter to another.</span></div>
<div class="viewcode-block" id="mpfit.tie"><a class="viewcode-back" href="../../agpy.html#agpy.mpfit_custom.mpfit.tie">[docs]</a>    <span class="k">def</span> <span class="nf">tie</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">ptied</span><span class="o">=</span><span class="bp">None</span><span class="p">):</span>
        <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">debug</span><span class="p">):</span> <span class="k">print</span> <span class="s">&#39;Entering tie...&#39;</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">ptied</span> <span class="o">==</span> <span class="bp">None</span><span class="p">):</span> <span class="k">return</span>
        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">ptied</span><span class="p">)):</span>
            <span class="k">if</span> <span class="n">ptied</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">==</span> <span class="s">&#39;&#39;</span><span class="p">:</span> <span class="k">continue</span>
            <span class="n">cmd</span> <span class="o">=</span> <span class="s">&#39;p[&#39;</span> <span class="o">+</span> <span class="nb">str</span><span class="p">(</span><span class="n">i</span><span class="p">)</span> <span class="o">+</span> <span class="s">&#39;] = &#39;</span> <span class="o">+</span> <span class="n">ptied</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
            <span class="k">exec</span><span class="p">(</span><span class="n">cmd</span><span class="p">)</span>
        <span class="k">return</span><span class="p">(</span><span class="n">p</span><span class="p">)</span>

    
    <span class="c">#    Original FORTRAN documentation</span>
    <span class="c">#    **********</span>
    <span class="c">#</span>
    <span class="c">#    subroutine covar</span>
    <span class="c">#</span>
    <span class="c">#    given an m by n matrix a, the problem is to determine</span>
    <span class="c">#    the covariance matrix corresponding to a, defined as</span>
    <span class="c">#</span>
    <span class="c">#                   t</span>
    <span class="c">#          inverse(a *a) .</span>
    <span class="c">#</span>
    <span class="c">#    this subroutine completes the solution of the problem</span>
    <span class="c">#    if it is provided with the necessary information from the</span>
    <span class="c">#    qr factorization, with column pivoting, of a. that is, if</span>
    <span class="c">#    a*p = q*r, where p is a permutation matrix, q has orthogonal</span>
    <span class="c">#    columns, and r is an upper triangular matrix with diagonal</span>
    <span class="c">#    elements of nonincreasing magnitude, then covar expects</span>
    <span class="c">#    the full upper triangle of r and the permutation matrix p.</span>
    <span class="c">#    the covariance matrix is then computed as</span>
    <span class="c">#</span>
    <span class="c">#                     t  t</span>
    <span class="c">#          p*inverse(r *r)*p  .</span>
    <span class="c">#</span>
    <span class="c">#    if a is nearly rank deficient, it may be desirable to compute</span>
    <span class="c">#    the covariance matrix corresponding to the linearly independent</span>
    <span class="c">#    columns of a. to define the numerical rank of a, covar uses</span>
    <span class="c">#    the tolerance tol. if l is the largest integer such that</span>
    <span class="c">#</span>
    <span class="c">#          abs(r(l,l)) .gt. tol*abs(r(1,1)) ,</span>
    <span class="c">#</span>
    <span class="c">#    then covar computes the covariance matrix corresponding to</span>
    <span class="c">#    the first l columns of r. for k greater than l, column</span>
    <span class="c">#    and row ipvt(k) of the covariance matrix are set to zero.</span>
    <span class="c">#</span>
    <span class="c">#    the subroutine statement is</span>
    <span class="c">#</span>
    <span class="c">#      subroutine covar(n,r,ldr,ipvt,tol,wa)</span>
    <span class="c">#</span>
    <span class="c">#    where</span>
    <span class="c">#</span>
    <span class="c">#      n is a positive integer input variable set to the order of r.</span>
    <span class="c">#</span>
    <span class="c">#      r is an n by n array. on input the full upper triangle must</span>
    <span class="c">#        contain the full upper triangle of the matrix r. on output</span>
    <span class="c">#        r contains the square symmetric covariance matrix.</span>
    <span class="c">#</span>
    <span class="c">#      ldr is a positive integer input variable not less than n</span>
    <span class="c">#        which specifies the leading dimension of the array r.</span>
    <span class="c">#</span>
    <span class="c">#      ipvt is an integer input array of length n which defines the</span>
    <span class="c">#        permutation matrix p such that a*p = q*r. column j of p</span>
    <span class="c">#        is column ipvt(j) of the identity matrix.</span>
    <span class="c">#</span>
    <span class="c">#      tol is a nonnegative input variable used to define the</span>
    <span class="c">#        numerical rank of a in the manner described above.</span>
    <span class="c">#</span>
    <span class="c">#      wa is a work array of length n.</span>
    <span class="c">#</span>
    <span class="c">#    subprograms called</span>
    <span class="c">#</span>
    <span class="c">#      fortran-supplied ... dabs</span>
    <span class="c">#</span>
    <span class="c">#    argonne national laboratory. minpack project. august 1980.</span>
    <span class="c">#    burton s. garbow, kenneth e. hillstrom, jorge j. more</span>
    <span class="c">#</span>
    <span class="c">#    **********</span>
    </div>
<div class="viewcode-block" id="mpfit.calc_covar"><a class="viewcode-back" href="../../agpy.html#agpy.mpfit_custom.mpfit.calc_covar">[docs]</a>    <span class="k">def</span> <span class="nf">calc_covar</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">rr</span><span class="p">,</span> <span class="n">ipvt</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">tol</span><span class="o">=</span><span class="mf">1.e-14</span><span class="p">):</span>

        <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">debug</span><span class="p">):</span> <span class="k">print</span> <span class="s">&#39;Entering calc_covar...&#39;</span>
        <span class="k">if</span> <span class="n">rank</span><span class="p">(</span><span class="n">rr</span><span class="p">)</span> <span class="o">!=</span> <span class="mi">2</span><span class="p">:</span>
            <span class="k">print</span> <span class="s">&#39;ERROR: r must be a two-dimensional matrix&#39;</span>
            <span class="k">return</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">)</span>
        <span class="n">s</span> <span class="o">=</span> <span class="n">shape</span><span class="p">(</span><span class="n">rr</span><span class="p">)</span>
        <span class="n">n</span> <span class="o">=</span> <span class="n">s</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
        <span class="k">if</span> <span class="n">s</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">!=</span> <span class="n">s</span><span class="p">[</span><span class="mi">1</span><span class="p">]:</span>
            <span class="k">print</span> <span class="s">&#39;ERROR: r must be a square matrix&#39;</span>
            <span class="k">return</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">)</span>

        <span class="k">if</span> <span class="p">(</span><span class="n">ipvt</span> <span class="o">==</span> <span class="bp">None</span><span class="p">):</span> <span class="n">ipvt</span> <span class="o">=</span> <span class="n">arange</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
        <span class="n">r</span> <span class="o">=</span> <span class="n">rr</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
        <span class="n">r</span><span class="o">.</span><span class="n">shape</span> <span class="o">=</span> <span class="p">[</span><span class="n">n</span><span class="p">,</span><span class="n">n</span><span class="p">]</span>

        <span class="c">## For the inverse of r in the full upper triangle of r</span>
        <span class="n">l</span> <span class="o">=</span> <span class="o">-</span><span class="mi">1</span>
        <span class="n">tolr</span> <span class="o">=</span> <span class="n">tol</span> <span class="o">*</span> <span class="nb">abs</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">])</span>
        <span class="k">for</span> <span class="n">k</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
            <span class="k">if</span> <span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">r</span><span class="p">[</span><span class="n">k</span><span class="p">,</span><span class="n">k</span><span class="p">])</span> <span class="o">&lt;=</span> <span class="n">tolr</span><span class="p">):</span> <span class="k">break</span>
            <span class="n">r</span><span class="p">[</span><span class="n">k</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.</span><span class="o">/</span><span class="n">r</span><span class="p">[</span><span class="n">k</span><span class="p">,</span><span class="n">k</span><span class="p">]</span>
            <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">k</span><span class="p">):</span>
                <span class="n">temp</span> <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="n">k</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">*</span> <span class="n">r</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">]</span>
                <span class="n">r</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.</span>
                <span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">-</span> <span class="n">temp</span><span class="o">*</span><span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span><span class="n">j</span><span class="p">]</span>
            <span class="n">l</span> <span class="o">=</span> <span class="n">k</span>

        <span class="c">## Form the full upper triangle of the inverse of (r transpose)*r</span>
        <span class="c">## in the full upper triangle of r</span>
        <span class="k">if</span> <span class="n">l</span> <span class="o">&gt;=</span> <span class="mi">0</span><span class="p">:</span>
            <span class="k">for</span> <span class="n">k</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">l</span><span class="o">+</span><span class="mi">1</span><span class="p">):</span>
                <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">k</span><span class="p">):</span>
                    <span class="n">temp</span> <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">]</span>
                    <span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">+</span> <span class="n">temp</span><span class="o">*</span><span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span><span class="n">k</span><span class="p">]</span>
                <span class="n">temp</span> <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="n">k</span><span class="p">,</span><span class="n">k</span><span class="p">]</span>
                <span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">k</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">=</span> <span class="n">temp</span> <span class="o">*</span> <span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">k</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span><span class="n">k</span><span class="p">]</span>

        <span class="c">## For the full lower triangle of the covariance matrix</span>
        <span class="c">## in the strict lower triangle or and in wa</span>
        <span class="n">wa</span> <span class="o">=</span> <span class="n">repeat</span><span class="p">([</span><span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">]],</span> <span class="n">n</span><span class="p">)</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
            <span class="n">jj</span> <span class="o">=</span> <span class="n">ipvt</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
            <span class="n">sing</span> <span class="o">=</span> <span class="n">j</span> <span class="o">&gt;</span> <span class="n">l</span>
            <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">):</span>
                <span class="k">if</span> <span class="n">sing</span><span class="p">:</span> <span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.</span>
                <span class="n">ii</span> <span class="o">=</span> <span class="n">ipvt</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
                <span class="k">if</span> <span class="n">ii</span> <span class="o">&gt;</span> <span class="n">jj</span><span class="p">:</span> <span class="n">r</span><span class="p">[</span><span class="n">ii</span><span class="p">,</span><span class="n">jj</span><span class="p">]</span> <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span>
                <span class="k">if</span> <span class="n">ii</span> <span class="o">&lt;</span> <span class="n">jj</span><span class="p">:</span> <span class="n">r</span><span class="p">[</span><span class="n">jj</span><span class="p">,</span><span class="n">ii</span><span class="p">]</span> <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span>
            <span class="n">wa</span><span class="p">[</span><span class="n">jj</span><span class="p">]</span> <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">j</span><span class="p">]</span>

        <span class="c">## Symmetrize the covariance matrix in r</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
            <span class="n">r</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">r</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="mi">0</span><span class="p">:</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span>
            <span class="n">r</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">wa</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>

        <span class="k">return</span><span class="p">(</span><span class="n">r</span><span class="p">)</span>
</div></div>
<div class="viewcode-block" id="machar"><a class="viewcode-back" href="../../agpy.html#agpy.mpfit_custom.machar">[docs]</a><span class="k">class</span> <span class="nc">machar</span><span class="p">:</span>
    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">double</span><span class="o">=</span><span class="mi">1</span><span class="p">):</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">double</span> <span class="o">==</span> <span class="mi">0</span><span class="p">):</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">machep</span> <span class="o">=</span> <span class="mf">1.19209e-007</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">maxnum</span> <span class="o">=</span> <span class="mf">3.40282e+038</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">minnum</span> <span class="o">=</span> <span class="mf">1.17549e-038</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">maxgam</span> <span class="o">=</span> <span class="mf">171.624376956302725</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">machep</span> <span class="o">=</span> <span class="mf">2.2204460e-016</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">maxnum</span> <span class="o">=</span> <span class="mf">1.7976931e+308</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">minnum</span> <span class="o">=</span> <span class="mf">2.2250739e-308</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">maxgam</span> <span class="o">=</span> <span class="mf">171.624376956302725</span>

        <span class="bp">self</span><span class="o">.</span><span class="n">maxlog</span> <span class="o">=</span> <span class="n">log</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">maxnum</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">minlog</span> <span class="o">=</span> <span class="n">log</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">minnum</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">rdwarf</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">minnum</span><span class="o">*</span><span class="mf">1.5</span><span class="p">)</span> <span class="o">*</span> <span class="mi">10</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">rgiant</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">maxnum</span><span class="p">)</span> <span class="o">*</span> <span class="mf">0.1</span>
</pre></div></div>

          </div>
        </div>
      </div>
        </div>
        <div class="sidebar">
          <h3>Table Of Contents</h3>
          <ul>
<li class="toctree-l1"><a class="reference internal" href="../../agpy.html">Adam Ginsburg&#8217;s Python Code (agpy)</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../image_tools.html">Image Tools</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../fft_tools.html">AG_fft_tools Package</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../plfit.html">plfit Package</a></li>
</ul>

          <h3 style="margin-top: 1.5em;">Search</h3>
          <form class="search" action="../../search.html" method="get">
            <input type="text" name="q" />
            <input type="submit" value="Go" />
            <input type="hidden" name="check_keywords" value="yes" />
            <input type="hidden" name="area" value="default" />
          </form>
          <p class="searchtip" style="font-size: 90%">
            Enter search terms or a module, class or function name.
          </p>
        </div>
        <div class="clearer"></div>
      </div>
    </div>

    <div class="footer">
      &copy; Copyright 2011, Adam Ginsburg.
      Created using <a href="http://sphinx.pocoo.org/">Sphinx</a> 1.2pre.
    <script type="text/javascript">

      var _gaq = _gaq || [];
      _gaq.push(['_setAccount', 'UA-6253248-2']);
      _gaq.push(['_trackPageview']);

      (function() {
        var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true;
        ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js';
        var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s);
      })();

    </script>
        
    </div>
  </body>
</html>